1 from pylab import *
   2 from scipy import *
   3 import time,calendar,os,pipes,struct,string
   4 import random
   5 import scipy.ndimage as ndi
   6 
   7 nk=3 #number of classes
   8 ndim=[100,100] #dimension of ASAR image in pixels
   9 
  10 fx0,fy0,fx2,fy2=[20.,100.,100.,10.] # freeboard coordinates from input file
  11 nf=11 #number of freeboard footprints
  12 fz_h=zeros(nf) #fz_h is the freeboard heights
  13 # produce a random array of freeboard heights
  14 for i in range(nf):
  15    fz_h[i]=random.random()*70
  16 
  17 # the x,y positions of the freeboard footprints
  18 if fx2==fx0:
  19    fx=zeros(nf)+fx0
  20    fy=linspace(0,ndim[1],nf)
  21 else:
  22    teta=(fy2-fy0)/(fx2-fx0)
  23    fx=linspace(fx0,fx2,nf)
  24    fy=(fx-fx0)*teta+fy0
  25 #produce a random matrix of ASAR image classes (in 3 classes)
  26 
  27 a=zeros((ndim[0],ndim[1]))
  28 for i in range(ndim[0]):
  29    for j in range(ndim[1]):
  30       a[i,j]=int(random.random()*(nk)+1)
  31 
  32 #tranform from float to int
  33 
  34 #A=ndimage.gaussian_filter(a,1)+0.5
  35 #B=A.astype(int)
  36 B=a.astype(int)
  37 
  38 #select the class of footprints from ASAR Image
  39 fz=zeros(nf)
  40 for i in range(nf):
  41    fz[i]=B[int(fx[i]-1),int(fy[i]-1)]
  42 
  43 #calculation of mean and standard deviation of freeboard heights in each class
  44 mean=zeros(nk)
  45 std_kl=zeros(nk)
  46 for kl in (1,2,3):
  47    s=(fz==kl)
  48    fh_kl=fz_h[s]
  49    mean[kl-1]=fh_kl.sum()/fh_kl.shape[0]
  50    std_kl[kl-1]=std(fh_kl)
  51 
  52 #plot mean and error bar
  53 kl=[1,2,3]
  54 bar(kl,mean,yerr=std_kl,ecolor='r',align='center')
  55 
  56 
  57 #calculate correlation coefficient
  58 #corr=xcorr(fz,fz_h,normed=True)[1][nf]
  59 
  60 figure()
  61 #plot ASAR image class together with freeboard heights
  62 x=arange(nf)
  63 bar(x,fz,align='center')
  64 ylim(0,6)
  65 ax=twinx()
  66 plot(x,fz_h,'r+-')
  67 xlim(-2,12)
  68 
  69 #plot ASAR image with freeboard flight
  70 figure()
  71 plot(fx,fy,'r+-',linewidth=2)
  72 xlim(0,100)
  73 ylim(0,100)
  74 colorbar()
  75 imshow(B,origin='lower',interpolation=None)
  76 hold(True)
  77 
  78 #figure()
  79 #imshow(a)
  80 show()