|
Size: 3116
Comment:
|
Size: 5618
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 14: | Line 14: |
| import string | from read_icesat import * from LinearAlgebra import * from Numeric import * |
| Line 16: | Line 18: |
| def fit_freeboard_ASAR(filename1,filename2,resolution): """filename1: ASAR data file, filename2: freeboard data file, resolution: data resolution""" |
def fit_freeboard_ASAR(filename1,filename2): """filename1: ASAR data file, filename2: freeboard data file creates new coordinate system defined by corners of ASAR image and selects freeboard values within ASAR image box returns an array containing normalized image coordinates and corresponding freeboard values: [x_coordinate, y_coordinate, freeboardheight(cm)]""" |
| Line 20: | Line 25: |
| lat1,lon1,lat2,lon2,lat3,lon3,lat4,lon4=read_asar_corners(filename1) ASAR_1=[lat1,lon1] ASAR_2=[lat2,lon2] ASAR_3=[lat3,lon3] ASAR_4=[lat4,lon4] ASAR_1p=mapll(ASAR_1[0],ASAR_1[1],sgn) #computing polarstereographic coordinates ASAR_2p=mapll(ASAR_2[0],ASAR_2[1],sgn) ASAR_3p=mapll(ASAR_3[0],ASAR_3[1],sgn) ASAR_4p=mapll(ASAR_4[0],ASAR_4[1],sgn) X=int(abs(ASAR_2p[0]-ASAR_1p[0])/resolution) #image size in pixel Y=int(abs(ASAR_4p[1]-ASAR_1p[1])/resolution) |
ASAR=array(read_asar_corners(filename1)) ASAR_p=zeros(8) for k in arange(0,7,2): #computing polarstereographic coordinates ASAR_p[k:k+2]=mapll(ASAR[k],ASAR[k+1],sgn) |
| Line 36: | Line 31: |
| 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]) | 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]) |
| Line 48: | Line 43: |
| # reading freeboard data lon=[] lat=[] fbh=[] datei = open (filename2, 'r') line=datei.readline() k=-1 while line!="": k=k+1 data=string.split(line) lon.append(float(data[0])) lat.append(abs(float(data[1]))) fbh.append(float(data[2])) line=datei.readline() polar=mapll(array(lat),array(lon),sgn) |
# reading freeboard data and computing geographic into polarstereographic coordinates ICESAT_p,fbh=read_icesat(filename2,sgn) |
| Line 68: | Line 49: |
| for x,y in zip(polar[0],polar[1]): | for x,y in zip(ICESAT_p[0],ICESAT_p[1]): |
| Line 72: | Line 53: |
| # cutting off non-corresponding data values | # cutting off non-corresponding data values m=-1 index_vec=[] for xn,yn in zip(x_neu,y_neu): m=m+1 if xn<=1. and xn >=0. and yn<=1. and yn >=0.: index_vec.append(m) x_bild=[] y_bild=[] fbh_bild=[] for i in index_vec: x_bild.append(x_neu[i]) y_bild.append(y_neu[i]) fbh_bild.append(fbh[i]) x_y_fbh=array([x_bild,y_bild,fbh_bild]) return x_y_fbh }}} {{{#!python import string from geo_polar import * from scipy import io def read_icesat(filename,sgn): data=io.read_array(filename) polar=mapll(data[:,1],data[:,0],sgn) fbh=data[:,2] return polar,fbh }}} Die benötigten Module polar_projection.py und read_asar.py sind auf der Seite der Arbeitsgruppe 0 [[AG0_ASAR_Einlesen]] zu finden. '''fbh_bildkoordinaten_test.py''' {{{#!python from polar_projection import * from read_asar import * from read_icesat import * from LinearAlgebra import * from Numeric import * def fit_freeboard_ASAR(filename1,filename2): """filename1: ASAR data file, filename2: freeboard data file creates new coordinate system defined by corners of ASAR image and selects freeboard values within ASAR image box returns an array containing normalized image coordinates and corresponding freeboard values: [x_coordinate, y_coordinate, freeboardheight(cm)]""" sgn=-1 #Antarctica ASAR=array(read_asar_corners(filename1)) ASAR_p=zeros(8) for k in arange(0,7,2): #computing polarstereographic coordinates ASAR_p[k:k+2]=mapll(ASAR[k],ASAR[k+1],sgn) # polarstereographic coordinate system 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]) # new coordinate system with normalized coordinates y10,x10,y11,x11,y12,x12,y13,x13=0,0,1,0,1,1,0,1 # calculating transformation matrix: P0=array([[x00, x01, x02, x03],[y00,y01,y02,y03],[1.0,1.0,1.0,1.0]]) P1=array([[x10, x11, x12, x13],[y10,y11,y12,y13],[1.0,1.0,1.0,1.0]]) Faktor1=dot(P1,transpose(P0)) Faktor2=inverse(dot(P0,transpose(P0))) A=dot(Faktor1,Faktor2) # Transformation matrix # reading freeboard data and computing geographic into polarstereographic coordinates ICESAT_p,fbh=read_icesat(filename2,sgn) # calculating new coordinates for freeboard data x_neu=[] y_neu=[] for x,y in zip(ICESAT_p[0],ICESAT_p[1]): x_neu.append(dot(array([A[0,0],A[0,1]]),array([x,y]))+A[0,2]) y_neu.append(dot(array([A[1,0],A[1,1]]),array([x,y]))+A[1,2]) # cutting off non-corresponding data values |
| Line 92: | Line 152: |
| filename1='ASA_IMP_1PNDPA20060617_043346_000000162048_00362_22460_2136.N1' filename2='LonLatFre_1706_6.xyz' ergebnis=fit_freeboard_ASAR(filename1,filename2) |
|
| Line 94: | Line 157: |
| Die benötigten Module polar_projection.py und read_asar.py sind auf der Seite zu finden. | {{attachment:schemabild.jpg}} |
Die Besprechung ueber das Vorgehen ihrer Aufgaben ergab:
- die Transformation in polarstereographische Koordinaten und der nachfolgende Uebergang zu den Bildpunkten von ASAR
Eckpunkte und Aufloseung des ASAR-Bildes variieren frei -> allgemeingueltiges Programm fuer den Uebergang der Output erfolgt als Vektor in der Form x, y, Freiboardhoehe; eventuell die Floatangabe
fbh_bildkoordinaten.py
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
Die benötigten Module polar_projection.py und read_asar.py sind auf der Seite der Arbeitsgruppe 0 AG0_ASAR_Einlesen zu finden.
fbh_bildkoordinaten_test.py
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
61
62 filename1='ASA_IMP_1PNDPA20060617_043346_000000162048_00362_22460_2136.N1'
63 filename2='LonLatFre_1706_6.xyz'
64 ergebnis=fit_freeboard_ASAR(filename1,filename2)
