1 from scipy import *
   2 from pyhdf.SD import *
   3 import os,os.path,string,pipes
   4 import struct#_hdfext
   5 from pylab import *
   6 
   7 def read_CHLO(filename):
   8         # Oeffne HDF4-Datei
   9         f=SD(filename)
  10         # Hole die Attribute (Metadata)
  11         attr = f.attributes()
  12         # Hole die in der Datei verfuegbaren Datensaetze
  13         dsets = f.datasets()
  14         # Hole die Daten aus der Datei in ein scipy.array
  15         data=array(f.select('l3m_data').get())
  16         # Siehe attr:
  17         #'Scaling Equation': ('Base**((Slope*l3m_data) + Intercept) = Parameter value\x00',
  18         Base=attr['Base']
  19         Slope=attr['Slope']
  20         Intercept=attr['Intercept']
  21         data[(data==65535).nonzero()]=nan
  22         # Skaliere die Daten, um die urspruenglichen Einheiten zu berechnen
  23         data=Base**((Slope*data) + Intercept)
  24         return data,attr
  25 
  26 def hdf_set(filename,z):
  27     hdfFile=filename
  28     f = SD(hdfFile,SDC.WRITE|SDC.CREATE)
  29     f.author = 'Lars Kaleschke'
  30     v2=f.create('z',SDC.FLOAT32,(z.shape[0],z.shape[1]))
  31     v2[:]=z.astype(float32)
  32     f.end()
  33     return
  34 
  35 def unzip_data(name1,data_dir):
  36         os.system('bunzip2 '+data_dir+name1)
  37         return
  38 
  39 def zip_data(name1,data_dir):
  40         os.system('bzip2 '+data_dir+name1)
  41         return
  42         
  43 def hdf_get(filename):
  44     f = SD(filename)
  45     dsets = f.datasets()
  46     d={}
  47     for dset in dsets.keys():
  48         ds=f.select(dset)
  49         d[dset]=array(ds.get(),float)
  50     return d
  51 
  52 
  53 def write_grd(data,lat,lon,filename,I,R):
  54 
  55     # Write to GMT/NetCDF .grd
  56 
  57     grdtable,grdfile=filename+'.tab',filename
  58     pipe=pipes.Template()
  59     cmd='xyz2grd -V  -bis '+I+' '+R+' -G'+grdfile
  60     print cmd
  61     pipe.append(cmd,'-.')
  62     f=pipe.open(grdtable, 'w')
  63     sy,sx=data.shape[0],data.shape[1]
  64     for yi in range(0,sy):
  65         for xi in range(0,sx):
  66             f.write(struct.pack("f",lon[yi,xi])+struct.pack("f",lat[yi,xi])+struct.pack("f",data[yi,xi]))
  67     f.close()
  68 
  69     return
  70 
  71 
  72 def gmt_map(data,outfile,title):
  73     Dlon=1.0/3
  74     Dlat=0.2
  75 
  76     X=88
  77     Y=82
  78 
  79     upper_lat=63.983
  80     lower_lat=upper_lat-Dlat*Y+Dlat
  81     western_lon=-15.25
  82     eastern_lon=western_lon+Dlon*X
  83 
  84     R=' -R%f'%western_lon+'/%f'%eastern_lon+'/%f'%lower_lat+'/%f'%upper_lat
  85     I=' -I0.3333333/0.2'
  86     J=' -Jm0.2 '
  87     A=' -Ba5g5/a5g5NESW '
  88     C='color.tab'
  89 
  90     lats=(upper_lat-indices((Y,X))[0,:,:].astype(float)*Dlat)
  91     lons=(indices((Y,X))[1,:,:].astype(float)*Dlon+western_lon)
  92 
  93     filename='data.grd'
  94     write_grd(data,lats,lons,filename,I,R)
  95 
  96 
  97     os.system('gmtset PAPER_MEDIA A4 PLOT_DEGREE_FORMAT ddd:mm:ss')
  98     os.system('makecpt -Cjet -T0/5/1  -Z -D  > '+C)
  99     os.system('grdimage '+filename+R+J+A+' -C'+C+'   -K > '+outfile)
 100     os.system('pscoast '+R+J+A+' -N1 -W1/0/0/0 -G200/200/200 -Di -K -O >> '+outfile)
 101 
 102     #Scale bar
 103     scl_tick='1/:mg/m\032:'
 104     print 'Add scale bar'
 105     os.system('psscale -D3.5/-0.6/2.0/0.1h -E -C'+C+' -B'+scl_tick+' -K -O >> '+outfile)
 106     os.system('echo "0.0 -0.03 14 0.0 0 0  '+title+'" | pstext -N -R0/1/0/1 -JX21/29  -O >> '+outfile)
 107     return
 108 
 109 
 110 for begin1 in range(1,365,8):
 111         datum1='2004%03d'%begin1
 112         end1=begin1+7
 113         datum2='2004%03d'%end1
 114         filename='A'+datum1+datum2+'.L3m_8D_CHLO_9.bz2'
 115         data_directory='/users/ifmlinux30a/ifmrs/u242023/SATBILD/data/'
 116         unzip_data(filename,data_directory)
 117         a=string.split((os.path.basename(data_directory+filename)),'.') #Aufspaltung des Dateinamenstrings in drei Teile
 118         b=a[0]+'.'+a[1]
 119         print b
 120 
 121         img,metadata=read_CHLO(data_directory+b) #Einlesen der Hdf-Daten
 122 
 123         ##ECOHAM3/4-Gebiet ausschneiden
 124         ##ECOHAM3/4 Latitude:47.5833 to 63.9833, Longitude:-15.25 to 14.0833
 125         ecolatstep=0.2
 126         ecolonstep=1./3.
 127         ecolat1=63.9833
 128         ecolat2=47.5833
 129         ecolon1=-15.25
 130         ecolon2=14.0833
 131 
 132         #Geographische Koordinaten des Satbildes:
 133         latstep=float(metadata['Latitude Step'])
 134         lonstep=float(metadata['Longitude Step'])
 135         Nlat=float(metadata['Northernmost Latitude'])
 136         Wlon=float(metadata['Westernmost Longitude'])
 137         Slat=float(metadata['Southernmost Latitude'])
 138         Olon=float(metadata['Easternmost Longitude'])
 139 
 140         #Berechnung der gedachten Modellraender bei Erweiterung auf (fast) gesamte Erde:
 141         modrand_N=ecolat1+floor((Nlat-ecolat1)/ecolatstep)*ecolatstep
 142         modrand_W=ecolon1+ceil((Wlon-ecolon1)/ecolonstep)*ecolonstep
 143         modrand_S=ecolat1+ceil((Slat-ecolat1)/ecolatstep)*ecolatstep
 144         modrand_O=ecolon1+floor((Olon-ecolon1)/ecolonstep)*ecolonstep
 145 
 146         #Berechnung der neuen Unterteilung des Satbildgitters fuer optimales Abschneiden zur Anpassung an Modellraender:
 147 
 148         zoom11=abs(latstep/(Nlat-modrand_N))   #Berechnung des zoom-Faktors, so dass am Nord- und Ostrand genau eine Pixelweite
 149         zoom12=abs(lonstep/(Wlon-modrand_W))   #abgeschnittten wird
 150         Sfaktor=int(round(abs(Slat-modrand_S)/(Nlat-modrand_N)))
 151         Ofaktor=int(round((Olon-modrand_O)/abs(Wlon-modrand_W)))
 152 
 153         img_zoom1=ndimage.zoom(img,(zoom11,zoom12),order=1)
 154 
 155         #Abschneiden der ueberstehenden Pixel im Satbild:
 156         s=img_zoom1.shape
 157         s1=s[0]-Sfaktor-1
 158         s2=s[1]-Ofaktor-1
 159         img_cut1=img_zoom1[1:s1,1:s2]
 160 
 161         #Anpassung des Satellitenbildes an Modellgitter:
 162         zoom21=(180./s[0])/0.2
 163         zoom22=(360./s[1])/(1./3.)
 164         img_zoom2=ndimage.zoom(img_cut1,(zoom21,zoom22),order=1)
 165 
 166         #Ausschneiden des Modellausschnitts aus Satbild:
 167         lat1=int(round((modrand_N-ecolat1)/ecolatstep))
 168         lat2=int(round((modrand_N-ecolat2)/ecolatstep))
 169         lon1=int(round(abs(modrand_W-ecolon1)/ecolonstep))
 170         lon2=int(round(abs(modrand_W-ecolon2)/ecolonstep))
 171         img_cut2_ord1=img_zoom2[lat1:lat2,lon1:lon2]
 172 
 173          
 174         #
 175         e=a[0]+'.'+a[1]+'_cut.hdf'
 176         hdf_set(data_directory+e,img_cut2_ord1) #Speichern des angepassten Satbildes als hdf-Datei
 177         #figure(1)  #Plotten
 178         #imshow(nan_to_num(img_cut2_ord1),vmin=0,vmax=5)
 179         #show()
 180         title='Chlorophyll '+a[0]
 181         outfile=a[0]+'.'+a[1]+'_cut.ps'
 182         gmt_map(img_cut2_ord1,data_directory+outfile,title)
 183 
 184         zip_data(b,data_directory) #komprimieren der gedownloadeten Dateien
 185         c=b+'.bz2'
 186         print c