1 from polar_projection import *
   2 from read_asar import *
   3 from read_icesat import *
   4 from LinearAlgebra import *
   5 from Numeric 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     # polarstereographic coordinate system
  20     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])
  21     # new coordinate system with normalized coordinates  
  22     y10,x10,y11,x11,y12,x12,y13,x13=0,0,1,0,1,1,0,1
  23 
  24     # calculating transformation matrix:
  25     P0=array([[x00, x01, x02, x03],[y00,y01,y02,y03],[1.0,1.0,1.0,1.0]])
  26     P1=array([[x10, x11, x12, x13],[y10,y11,y12,y13],[1.0,1.0,1.0,1.0]])
  27 
  28     Faktor1=dot(P1,transpose(P0))
  29     Faktor2=inverse(dot(P0,transpose(P0)))
  30     A=dot(Faktor1,Faktor2)  # Transformation matrix 
  31 
  32     # reading freeboard data and computing geographic into polarstereographic coordinates
  33     ICESAT_p,fbh=read_icesat(filename2,sgn)
  34 
  35     # calculating new coordinates for freeboard data
  36     x_neu=[]
  37     y_neu=[]
  38     for x,y in zip(ICESAT_p[0],ICESAT_p[1]):
  39         x_neu.append(dot(array([A[0,0],A[0,1]]),array([x,y]))+A[0,2])
  40         y_neu.append(dot(array([A[1,0],A[1,1]]),array([x,y]))+A[1,2])
  41 
  42     # cutting off non-corresponding data values  
  43     m=-1
  44     index_vec=[]
  45     for xn,yn in zip(x_neu,y_neu):
  46         m=m+1
  47         if xn<=1. and xn >=0. and yn<=1. and yn >=0.:
  48             index_vec.append(m)
  49 
  50     x_bild=[]
  51     y_bild=[] 
  52     fbh_bild=[]
  53     for i in index_vec:
  54         x_bild.append(x_neu[i])
  55         y_bild.append(y_neu[i])
  56         fbh_bild.append(fbh[i])
  57 
  58     x_y_fbh=array([x_bild,y_bild,fbh_bild])
  59 
  60     return x_y_fbh 

   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)