SiaChlorophyllProjekt

=== Die einzelnen Teilergebnisse der Gruppen sind hier zusammen gefügt: ===
Erklärungen zu den einzelnen Programmteilen sind unter ChlorophyllArbeitsGruppe1 und ChlorophyllArbeitsGruppe2 zu finden. Die Funktionen und Programmteile, die die gmt-Darstellung ermöglichen wurden von Lars Kaleschke erarbeitet (siehe auch ChlorophyllProgrammcodes ).


{{{#!python
from scipy import *
from pyhdf.SD import *
import os,os.path,string,pipes
import struct#_hdfext
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

def unzip_data(name1,data_dir):
        os.system('bunzip2 '+data_dir+name1)
        return

def zip_data(name1,data_dir):
        os.system('bzip2 '+data_dir+name1)
        return
	
def hdf_get(filename):
    f = SD(filename)
    dsets = f.datasets()
    d={}
    for dset in dsets.keys():
        ds=f.select(dset)
        d[dset]=array(ds.get(),float)
    return d


def write_grd(data,lat,lon,filename,I,R):

    # Write to GMT/NetCDF .grd

    grdtable,grdfile=filename+'.tab',filename
    pipe=pipes.Template()
    cmd='xyz2grd -V  -bis '+I+' '+R+' -G'+grdfile
    print cmd
    pipe.append(cmd,'-.')
    f=pipe.open(grdtable, 'w')
    sy,sx=data.shape[0],data.shape[1]
    for yi in range(0,sy):
        for xi in range(0,sx):
            f.write(struct.pack("f",lon[yi,xi])+struct.pack("f",lat[yi,xi])+struct.pack("f",data[yi,xi]))
    f.close()

    return


def gmt_map(data,outfile,title):
    Dlon=1.0/3
    Dlat=0.2

    X=88
    Y=82

    upper_lat=63.983
    lower_lat=upper_lat-Dlat*Y+Dlat
    western_lon=-15.25
    eastern_lon=western_lon+Dlon*X

    R=' -R%f'%western_lon+'/%f'%eastern_lon+'/%f'%lower_lat+'/%f'%upper_lat
    I=' -I0.3333333/0.2'
    J=' -Jm0.2 '
    A=' -Ba5g5/a5g5NESW '
    C='color.tab'

    lats=(upper_lat-indices((Y,X))[0,:,:].astype(float)*Dlat)
    lons=(indices((Y,X))[1,:,:].astype(float)*Dlon+western_lon)

    filename='data.grd'
    write_grd(data,lats,lons,filename,I,R)


    os.system('gmtset PAPER_MEDIA A4 PLOT_DEGREE_FORMAT ddd:mm:ss')
    os.system('makecpt -Cjet -T0/5/1  -Z -D  > '+C)
    os.system('grdimage '+filename+R+J+A+' -C'+C+'   -K > '+outfile)
    os.system('pscoast '+R+J+A+' -N1 -W1/0/0/0 -G200/200/200 -Di -K -O >> '+outfile)

    #Scale bar
    scl_tick='1/:mg/m\032:'
    print 'Add scale bar'
    os.system('psscale -D3.5/-0.6/2.0/0.1h -E -C'+C+' -B'+scl_tick+' -K -O >> '+outfile)
    os.system('echo "0.0 -0.03 14 0.0 0 0  '+title+'" | pstext -N -R0/1/0/1 -JX21/29  -O >> '+outfile)
    return


for begin1 in range(1,365,8):
        datum1='2004%03d'%begin1
        end1=begin1+7
        datum2='2004%03d'%end1
        filename='A'+datum1+datum2+'.L3m_8D_CHLO_9.bz2'
        data_directory='/users/ifmlinux30a/ifmrs/u242023/SATBILD/data/'
        unzip_data(filename,data_directory)
        a=string.split((os.path.basename(data_directory+filename)),'.') #Aufspaltung des Dateinamenstrings in drei Teile
        b=a[0]+'.'+a[1]
        print b

        img,metadata=read_CHLO(data_directory+b) #Einlesen der Hdf-Daten

	##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

	#Geographische Koordinaten des Satbildes:
	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'])

	#Berechnung der gedachten Modellraender bei Erweiterung auf (fast) gesamte Erde:
	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 zur Anpassung an Modellraender:

	zoom11=abs(latstep/(Nlat-modrand_N))   #Berechnung des zoom-Faktors, so dass am Nord- und Ostrand genau eine Pixelweite
	zoom12=abs(lonstep/(Wlon-modrand_W))   #abgeschnittten wird
	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),order=1)

	#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),order=1)

	#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_ord1=img_zoom2[lat1:lat2,lon1:lon2]

         
	#
        e=a[0]+'.'+a[1]+'_cut.hdf'
        hdf_set(data_directory+e,img_cut2_ord1) #Speichern des angepassten Satbildes als hdf-Datei
        #figure(1)  #Plotten
        #imshow(nan_to_num(img_cut2_ord1),vmin=0,vmax=5)
        #show()
	title='Chlorophyll '+a[0]
	outfile=a[0]+'.'+a[1]+'_cut.ps'
	gmt_map(img_cut2_ord1,data_directory+outfile,title)

        zip_data(b,data_directory) #komprimieren der gedownloadeten Dateien
        c=b+'.bz2'
        print c
}}}