1 from geo_polar import *
   2 from read_asar import *
   3 import string
   4 
   5 def fit_freeboard_ASAR(filename1,filename2,resolution):
   6     """filename1: ASAR data file, filename2: freeboard data file, resolution: data resolution"""
   7     
   8     sgn=-1  #Antarctica
   9     lat1,lon1,lat2,lon2,lat3,lon3,lat4,lon4=read_asar_corners(filename1)
  10     
  11     ASAR_1=[lat1,lon1]      
  12     ASAR_2=[lat2,lon2]      
  13     ASAR_3=[lat3,lon3]      
  14     ASAR_4=[lat4,lon4]     
  15 
  16     ASAR_1p=mapll(ASAR_1[0],ASAR_1[1],sgn)  #computing polarstereographic coordinates
  17     ASAR_2p=mapll(ASAR_2[0],ASAR_2[1],sgn)
  18     ASAR_3p=mapll(ASAR_3[0],ASAR_3[1],sgn)
  19     ASAR_4p=mapll(ASAR_4[0],ASAR_4[1],sgn)
  20 
  21     X=int(abs(ASAR_2p[0]-ASAR_1p[0])/resolution)  #image size in pixel
  22     Y=int(abs(ASAR_4p[1]-ASAR_1p[1])/resolution)
  23 
  24     # polarstereographic coordinate system
  25     y00,x00,y01,x01,y02,x02,y03,x03=int(ASAR_1p[1]),int(ASAR_1p[0]),int(ASAR_4p[1]),int(ASAR_4p[0]),int(ASAR_3p[1]),int(ASAR_3p[0]),int(ASAR_2p[1]),int(ASAR_2p[0])
  26     # new coordinate system with normalized coordinates  
  27     y10,x10,y11,x11,y12,x12,y13,x13=0,0,1,0,1,1,0,1
  28 
  29     # calculating transformation matrix:
  30     P0=array([[x00, x01, x02, x03],[y00,y01,y02,y03],[1.0,1.0,1.0,1.0]])
  31     P1=array([[x10, x11, x12, x13],[y10,y11,y12,y13],[1.0,1.0,1.0,1.0]])
  32 
  33     Faktor1=dot(P1,transpose(P0))
  34     Faktor2=inverse(dot(P0,transpose(P0)))
  35     A=dot(Faktor1,Faktor2)  # Transformation matrix 
  36 
  37     # reading freeboard data
  38     lon=[]
  39     lat=[]
  40     fbh=[]
  41     datei = open (filename2, 'r')
  42     line=datei.readline()
  43     k=-1
  44     while line!="":
  45         k=k+1
  46         data=string.split(line)
  47         lon.append(float(data[0]))
  48         lat.append(abs(float(data[1])))
  49         fbh.append(float(data[2]))
  50         line=datei.readline()
  51 
  52     polar=mapll(array(lat),array(lon),sgn)
  53 
  54     # calculating new coordinates for freeboard data
  55     x_neu=[]
  56     y_neu=[]
  57     for x,y in zip(polar[0],polar[1]):
  58         x_neu.append(dot(array([A[0,0],A[0,1]]),array([x,y]))+A[0,2])
  59         y_neu.append(dot(array([A[1,0],A[1,1]]),array([x,y]))+A[1,2])
  60 
  61     # cutting off non-corresponding data values 
  62     m=-1
  63     index_vec=[]
  64     for xn,yn in zip(x_neu,y_neu):
  65         m=m+1
  66         if xn<=1. and xn >=0. and yn<=1. and yn >=0.:
  67             index_vec.append(m)
  68 
  69     x_bild=[]
  70     y_bild=[] 
  71     fbh_bild=[]
  72     for i in index_vec:
  73         x_bild.append(x_neu[i])
  74         y_bild.append(y_neu[i])
  75         fbh_bild.append(fbh[i])
  76 
  77     x_y_fbh=array([x_bild,y_bild,fbh_bild])
  78 
  79     return x_y_fbh