Differences between revisions 10 and 12 (spanning 2 versions)
Revision 10 as of 2008-07-10 09:01:46
Size: 6629
Editor: NinaMaass
Comment:
Revision 12 as of 2008-07-10 10:06:32
Size: 3409
Editor: NinaMaass
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: ASAR image data resolution
       creates new coordinate system defined by corners of ASAR image and selects freeboard values within
       ASAR image box and returns an array containing normalized image coordinates and corresponding
       freeboard values:
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:
Line 24: 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 40: 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 52: 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 72: Line 49:
    for x,y in zip(polar[0],polar[1]):     for x,y in zip(ICESAT_p[0],ICESAT_p[1]):
Line 76: Line 53:
    # cutting off non-corresponding data values     # cutting off non-corresponding data values  
Line 95: Line 72:
}}}
Line 96: Line 74:
{{{#!python
import string
from polar_projection 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
Line 102: Line 90:
Zum Testen hängt man an das obige Programm folgende Zeilen an:
Line 103: Line 92:
from polar_projection import *
from read_asar import *
import string

def fit_freeboard_ASAR(filename1,filename2,resolution):
    """filename1: ASAR data file, filename2: freeboard data file, resolution: ASAR image data resolution
       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
    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)

    # polarstereographic coordinate system
    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])
    # 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
    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)

    # calculating new coordinates for freeboard data
    x_neu=[]
    y_neu=[]
    for x,y in zip(polar[0],polar[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
    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
Line 189: Line 94:
resolution=0.025
ergebnis=fit_freeboard_ASAR(filename1,filename2,resolution)
ergebnis=fit_freeboard_ASAR(filename1,filename2)

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 

   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

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

Zum Testen hängt man an das obige Programm folgende Zeilen an:

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

schemabild.jpg

LehreWiki: \AG1_Freibord (last edited 2008-07-11 11:19:34 by NinaMaass)