SiaChlorophyllProjekt

== Arbeitsgruppe 2 : Nina Maaß, Bente Tiedje ==

'''''Ziel/Aufgabe:''''' Die heruntergeladenen Satellitendaten sollen so bearbeitet werden, dass sie mit der ECOHAM-Modellgitterstruktur (88x82 Boxen) übereinstimmen. Anschließend sollen sie auf das Modellgebiet (Nordseeregion) verkleinert werden.

==== Einfaches Einlesen, Ausschneiden und Darstellen ====
Mit diesem Programm werden die Daten (A20040012004366.L3m_YR_CHLO_4) eingelesen. Durch einfache Umrechnung der Zeilen und Spalten der "Datenmatrix" in Breiten- und Längengrade werden die am nächsten liegenden geographischen Eckpunkte des ECOHAM-Modellgebiet ermittelt. Anschließend wird das ECOHAM-Modellgebiet ausgeschnitten und dargestellt (siehe Abbildung 1). Da Rundungen vorgenommen werden müssen, wird nicht das exakte Modellgebiet ausgeschnitten. Eine weitaus genauere Berechnung des Modellgebietes liefert erst das nachfolgende Programm. 



{{{#!python
from scipy import * 
Abbildung 4
import gdal 
import struct
from pylab import *

filename='/pf/u/u243066/satdat/chlorophyll/A20040012004366.L3m_YR_CHLO_4' 
f=gdal.Open(filename)         
metadata=f.GetMetadata()       
a=f.GetRasterBand(1)          
scanline = a.ReadRaster( 0, 0, a.XSize, 1, a.XSize, 1, gdal.GDT_Float32 )

img=zeros((a.YSize,a.XSize),float)
for yi in range(a.YSize-1):
        scanline=struct.unpack("f"*a.XSize,a.ReadRaster(0,yi,a.XSize,1,a.XSize,1,gdal.GDT_Float32))
        img[yi,:]=((array(scanline).astype(float)).astype(float))

#img enthält Daten für gesamte Erde:90°N-90°S, 180°W-180°O

#damit die Landfläche dunkelblau erscheint, also in der "Farbe" von Null
img[(img==65535).nonzero()]=0 

##ECOHAM3/4-Gebiet ausschneiden
##ECOHAM3/4 Latitude:47.5833 to 63.9833, Longitude:-15.25 to 14.0833

ecolat1=63.9833
ecolat2=47.5833
ecolon1=-15.25
ecolon2=14.0833                               
latstep=float(metadata['Latitude Step'])     
lonstep=float(metadata['Longitude Step'])    
Nlat=float(metadata['Northernmost Latitude'])
Wlon=float(metadata['Westernmost Longitude'])

lat1=floor((Nlat-ecolat1)/latstep)#obere Breitengradbegrenzung des Gebietes
lat2=ceil((Nlat-ecolat2)/latstep) #untere Breitengradbegrenzung des Gebietes
lon1=floor((abs(Wlon-ecolon1))/lonstep)#linke Längengradbegrenzung des Gebietes
lon2=ceil((abs(Wlon-ecolon2))/lonstep)#rechte Längengradbegrenzung des Gebietes

#Ausschneiden des Gebietes
img_eco34=img[lat1:lat2,lon1:lon2]

#Skalierung der Daten (siehe key 'Scaling Equation' in metadata):
img_eco34_chlo=float(metadata['Base'])**((float(metadata['Slope'])*img_eco34)
+float(metadata['Intercept']))
#Plotten
figure(1)
imshow(img_eco34_chlo,vmin=0,vmax=5)
title('MODIS 2004 Chlorophyll 4km')
colorbar(shrink=0.6)
show()
}}}
Abbildung 1: Ergebnis-Plot des Programmcodes

{{attachment:A20040012004366.L3m_YR_CHLO_4.png}} 

==== Anpassen an Modellgitterstruktur (und Ausschneiden) ====

Um den nachfolgenden Programmcode besser nachvollziehen zu können, wird er jetzt kurz erklärt. Die Abbildungen 2 und 3 veranschaulichen verwendete Konstantenbezeichnungen.

 * Zunächst wird das Modellgitter (Gitterweite N-S: 0,2° und O-W: 0,333°) auf die ganze Erde ausgeweitet. Anschließend wird berechnet wie groß die Differenzen zwischen den Ausmaßen der Satellitendaten (90°N bis 90°S, 180°W bis 180°O) und den Ausmaßen des erweiterten Modellgebietes sind. Das Verhältnis der Differenzen vom linken zum rechten Rand beträgt 1:3 (vergleiche Ofaktor) und vom oberen zum unteren Rand 1:11 (vergleiche Sfaktor).
 * Die Satellitendaten werden nun mit Hilfe der zoom-Funktion so unterteilt (interpoliert), dass eine Schrittweite in N-S-Richtung der berechneteten Differenz des oberen Randes und eine Schrittweite in O-W-Richtung der berechneten Differenz des linken Randes entspricht (siehe Zeile 66).
 * Von den Satellitendaten werden die erste Zeile und erste Spalte, sowie die 11 letzten Zeilen und die 3 letzten Spalten "abgeschnitten" (Zeile 72). 
 * Auf diese übrig gebliebenen Daten (img_cut1) wird die zoom-Funktion so angewandt, dass die Schrittweiten den Gitterweiten des Modells entsprechen.
 * Die Eckpunkte des Modellgebietes (lat1, lat2, lon1, lon2) werden berechnet und das Modellgebiet wird ausgeschnitten (Zeile 79 bis 84). Die ausgeschnittenen Daten (img_cut2) haben nun genausoviele Zeilen und Spalten wie das Modell Boxen in lateraler und longitudinaler Richtung hat, nämlich 82x88.
 * img_cut2 wird als hdf-Datei gespeichert und zur Veranschaulichung graphisch dargestellt (siehe Abbildung 4).

{{attachment:Modis_2004_YR_4320x8640.png}} Abbildung 2

{{attachment:Modellgitter_82x88.jpg}} Abbildung 3

{{{#!python
from scipy import *
import struct
from pyhdf.SD import *
from pylab import *

def read_CHLO(filename):
        # Oeffne HDF4-Datei
        f=SD(filename)
        # Hole die Attribute (Metadata)
        attr = f.attributes()
        # Hole die in der Datei verfuegbaren Datensaetze
        dsets = f.datasets()
        # Hole die Daten aus der Datei in ein scipy.array
        data=array(f.select('l3m_data').get())
        # Siehe attr:
        #'Scaling Equation': ('Base**((Slope*l3m_data) + Intercept) = Parameter value\x00',
        Base=attr['Base']
        Slope=attr['Slope']
        Intercept=attr['Intercept']
        data[(data==65535).nonzero()]=nan
        # Skaliere die Daten, um die urspruenglichen Einheiten zu berechnen
        data=Base**((Slope*data) + Intercept)
        return data,attr

def hdf_set(filename,z):
    hdfFile=filename
    f = SD(hdfFile,SDC.WRITE|SDC.CREATE)
    f.author = 'Lars Kaleschke'
    v2=f.create('z',SDC.FLOAT32,(z.shape[0],z.shape[1]))
    v2[:]=z.astype(float32)
    f.end()
    return

filename='/pf/u/u243066/satdat/chlorophyll/A20040012004366.L3m_YR_CHLO_4'
img,metadata=read_CHLO(filename)

##ECOHAM3/4-Gebiet ausschneiden
##ECOHAM3/4 Latitude:47.5833 to 63.9833, Longitude:-15.25 to 14.0833

ecolatstep=0.2
ecolonstep=1./3.
ecolat1=63.9833
ecolat2=47.5833
ecolon1=-15.25
ecolon2=14.0833
latstep=float(metadata['Latitude Step'])
lonstep=float(metadata['Longitude Step'])
Nlat=float(metadata['Northernmost Latitude'])
Wlon=float(metadata['Westernmost Longitude'])
Slat=float(metadata['Southernmost Latitude'])
Olon=float(metadata['Easternmost Longitude'])

#Anpassung des Satellitenbildrandes an Modellrand:
modrand_N=ecolat1+floor((Nlat-ecolat1)/ecolatstep)*ecolatstep
modrand_W=ecolon1+ceil((Wlon-ecolon1)/ecolonstep)*ecolonstep
modrand_S=ecolat1+ceil((Slat-ecolat1)/ecolatstep)*ecolatstep
modrand_O=ecolon1+floor((Olon-ecolon1)/ecolonstep)*ecolonstep

#Berechnung der neuen Unterteilung des Satbildgitters fuer optimales Abschneiden:

zoom11=abs(latstep/(Nlat-modrand_N))   #Berechnung des zoom-Faktors
zoom12=abs(lonstep/(Wlon-modrand_W))
Sfaktor=int(round(abs(Slat-modrand_S)/(Nlat-modrand_N)))
Ofaktor=int(round((Olon-modrand_O)/abs(Wlon-modrand_W)))

img_zoom1=ndimage.zoom(img,(zoom11,zoom12)) #Default-Einstellung für Grad der zoom-Funktion: order=3

#Abschneiden der ueberstehenden Pixel im Satbild:
s=img_zoom1.shape
s1=s[0]-Sfaktor-1
s2=s[1]-Ofaktor-1
img_cut1=img_zoom1[1:s1,1:s2]

#Anpassung des Satellitenbildes an Modellgitter:
zoom21=(180./s[0])/0.2
zoom22=(360./s[1])/(1./3.)
img_zoom2=ndimage.zoom(img_cut1,(zoom21,zoom22)) #Default-Einstellung für Grad der zoom-Funktion: order=3

#Ausschneiden des Modellausschnitts aus Satbild:
lat1=int(round((modrand_N-ecolat1)/ecolatstep))
lat2=int(round((modrand_N-ecolat2)/ecolatstep))
lon1=int(round(abs(modrand_W-ecolon1)/ecolonstep))
lon2=int(round(abs(modrand_W-ecolon2)/ecolonstep))
img_cut2=img_zoom2[lat1:lat2,lon1:lon2]

#Speichern von img_cut2 als hdf-Datei:
filename='Satimg_cut_Eco34.hdf'
hdf_set(filename,img_cut2)

#Plotten von img_cut2:
figure(3)
imshow(nan_to_num(img_cut2),vmin=0,vmax=5)
show()
 
}}}

Abbildung 4: Ergebnis-Plot des Programmcodes
 
{{attachment:satbild_auf_modgit_82x88.jpg}}