1 from polar_projection import *
   2 from read_asar import *
   3 from read_icesat import *
   4 from coord_transform import *
   5 from scipy import *
   6 
   7 def fit_freeboard_ASAR(filename1,filename2):
   8     """filename1: ASAR data file, filename2: freeboard data file
   9        creates new coordinate system defined by corners of ASAR image and selects freeboard values within ASAR image box
  10        returns an array containing normalized image coordinates and corresponding freeboard values:
  11        [x_coordinate, y_coordinate, freeboardheight(cm)]"""
  12     
  13     sgn=-1  #Antarctica
  14     ASAR=array(read_asar_corners(filename1))
  15     ASAR_p=zeros(8)
  16     for k in arange(0,7,2):       #computing polarstereographic coordinates
  17         ASAR_p[k:k+2]=mapll(ASAR[k],ASAR[k+1],sgn)
  18 
  19     A=coord_transformation(ASAR_p)
  20 
  21     # reading freeboard data and computing geographic into polarstereographic coordinates
  22     ICESAT_p,fbh=read_icesat(filename2,sgn)    #fbh are measured freeboard heights in cm
  23 
  24     # calculating new coordinates for freeboard data
  25     x_neu=[]
  26     y_neu=[]
  27     for x,y in zip(ICESAT_p[0],ICESAT_p[1]):
  28         x_neu.append(dot(array([A[0,0],A[0,1]]),array([x,y]))+A[0,2])
  29         y_neu.append(dot(array([A[1,0],A[1,1]]),array([x,y]))+A[1,2])
  30    
  31     # cutting off non-corresponding data values  
  32     m=-1
  33     index_vec=[]
  34     for xn,yn in zip(x_neu,y_neu):
  35         m=m+1
  36         if xn<=1. and xn >=0. and yn<=1. and yn >=0.:
  37             index_vec.append(m)
  38 
  39     x_bild=[]
  40     y_bild=[] 
  41     fbh_bild=[]
  42     for i in index_vec:
  43         x_bild.append(x_neu[i])
  44         y_bild.append(y_neu[i])
  45         fbh_bild.append(fbh[i])
  46 
  47     x_y_fbh=array([x_bild,y_bild,fbh_bild])
  48 
  49     return x_y_fbh

   1 from Numeric import *
   2 import LinearAlgebra as la
   3 
   4 def coord_transformation(ASAR_p):
   5     # polarstereographic coordinate system
   6     y00,x00,y01,x01,y02,x02,y03,x03=int(ASAR_p[1]),int(ASAR_p[0]),int(ASAR_p[7]),int(ASAR_p[6]),int(ASAR_p[5]),int(ASAR_p[4]),int(ASAR_p[3]),int(ASAR_p[2])
   7     # new coordinate system with normalized coordinates  
   8     y10,x10,y11,x11,y12,x12,y13,x13=0,0,1,0,1,1,0,1
   9 
  10     # calculating transformation matrix:
  11     P0=array([[x00, x01, x02, x03],[y00,y01,y02,y03],[1.0,1.0,1.0,1.0]])
  12     P1=array([[x10, x11, x12, x13],[y10,y11,y12,y13],[1.0,1.0,1.0,1.0]])
  13 
  14     Faktor1=dot(P1,transpose(P0))
  15     Faktor2=la.inverse(dot(P0,transpose(P0)))
  16     A=dot(Faktor1,Faktor2)  # Transformation matrix
  17     return A

   1 import string
   2 from polar_projection import *
   3 from scipy import io
   4 
   5 def read_icesat(filename,sgn):
   6     data=io.read_array(filename)
   7     polar=mapll(data[:,1],data[:,0],sgn)
   8     fbh=data[:,2]
   9     return polar,fbh

   1 filename1='ASA_IMP_1PNDPA20060617_043346_000000162048_00362_22460_2136.N1'
   2 filename2='LonLatFre_1706_6.xyz'
   3 ergebnis=fit_freeboard_ASAR(filename1,filename2)