SiaChlorophyllProjekt

== Fuellen der Luecken ==
{{attachment:filling_gaps.png}}

{{{#!python
from scipy import *
from pylab import *
from scipy.ndimage import *

def interpolate_gaps(Z):
	ind_nan=isnan(Z) # Index of NaNs
	sy,sx=Z.shape[0],z.shape[1]
	Z0=Z.copy()
	ZY=Z.copy()
	ZX=Z.copy()
	
	xi=arange(sx)
	for y in range(sy):
		L=Z[y,:]
		ind=isfinite(L)
		X=xi[ind]
		Y=L[ind]
		tck=interpolate.splrep(X,Y,s=0.0) # Calculate coefficients at the supporting points
		spl_int=interpolate.splev(xi,tck)
		ZY[y,:]=spl_int
		
	Z0[ind_nan]=ZY[ind_nan]
	return Z0


z=rand(60,60)*10
z=gaussian_filter(z,sigma=3)
close('all')
figure(1)
z0=z.copy()# Original image
z1=z.copy()# Original with missing values (set to zero)
#z2 temporary dilated image
z3=z.copy()# Reconstructed image with filled values

subplot(221)
imshow(z,vmin=4,vmax=6,interpolation='nearest')
colorbar()
title('Original')

# Generate missing elements set to zero
for i in range(30):
	s=int(rand(1)*10)
	ry=int(rand(1)*z.shape[0])
	rx=int(rand(1)*z.shape[1])	
	z1[ry:ry+s,rx:rx+s]=nan

ind_zero=(z1==0).nonzero() # Index of elements that are zero

z2=z1.copy()
z3=interpolate_gaps(z2)


subplot(222)
imshow(z1,vmin=4,vmax=6,interpolation='nearest')
colorbar()
title('Missing data')

subplot(223)
imshow(z3,vmin=4,vmax=6,interpolation='nearest')
colorbar()
title('Filled gaps')

subplot(224)
imshow(z0-z3,interpolation='nearest')
colorbar()
title('Difference')
savefig('filling_gaps.png',dpi=50)
show()
}}}


== Erzeugung und Speicherung einer Testdatei für die gmt-Darstellung ==
Eine Erklärung zu diesem Programm findet man unter ChlorophyllArbeitsGruppe2 .

{{{#!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/Modis_2004_YR_chlo_9'
img,metadata=read_CHLO(filename)

##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)
lat2=ceil((Nlat-ecolat2)/latstep)
lon1=floor((abs(Wlon-ecolon1))/lonstep)
lon2=ceil((abs(Wlon-ecolon2))/lonstep)
img_eco34=img[lat1:lat2,lon1:lon2]

#Anpassung des Satellitenbildausschnittes auf Modellgitter:
img_eco34_modgit=ndimage.zoom(img_eco34,(82./198.,88./353.))

#Speichern des angepassten Satbildes als hdf-Datei:
 
filename='testbild82x88.hdf'
hdf_set(filename,img_eco34_modgit)
}}}

== gmt-Darstellung ==

{{{#!python
from scipy import *
import struct,pipes,os
from pyhdf.SD import *

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



#TEST
filename='testbild82x88.hdf'
d=hdf_get(filename)

data=d['z']
outfile='map.ps'
title='Chlorophyll am Tag X'
gmt_map(data,outfile,title)
os.system('gv '+outfile)
}}}

=== GMT Karte ===
{{attachment:gmt_map.png}}