1 #!/bin/tcsh
   2 
   3 from pyhdf import SD
   4 from pylab import *
   5 from mpl_toolkits.basemap import Basemap
   6 
   7 fid=SD.SD('AIRS.2003.03.30.L3.RetStd001.v5.0.14.0.G07215021620.hdf')
   8 T=fid.select('Temperature_A')
   9 Temp=T.get()
  10 
  11 
  12 #Drucklevel
  13 Lv=[1000.0,925.0,850.0,700.0,600.0,500.0,400.0,300.0,250.0,200.0,150.0,100.0,
  14   70.0, 50.0,30.0, 20.0, 15.0, 10.0,  7.0,5.0, 3.0,2.0,1.5,1.0]
  15 
  16 
  17 # Fehlende Profile finden - Wert 1 
  18 testnan=Temp[0,:,:]
  19 testnan=testnan<0
  20 
  21 # Temperaturprofile
  22 figure(1)
  23 subplot(1,2,1)
  24 plot(Temp[:,16,165],Lv)
  25 ylim((1000.0,1.0))
  26 title('over ice, lat=-15,lon=74')
  27 show()
  28 subplot(1,2,2)
  29 plot(Temp[:,16,192],Lv)
  30 ylim((1000.0,1.0))
  31 title('over open water, lat=12,lon=74')
  32 show()
  33 savefig('Soundings.png',dpi=75)
  34 
  35 # Kartendarstellung der vorhandenen Profile
  36 figure(2)
  37 m = Basemap(width=2400000,height=1600000,projection='stere',lat_ts=77.0,lon_0=15.0,lat_0=77.0,resolution='l')
  38 m.drawcoastlines()
  39 
  40 for i in range(0,30):
  41     for j in range(160,200):
  42         if testnan[i,j]==0:
  43             x,y=m(j-180.,90.-i)
  44             m.plot(x,y,'r.')
  45 
  46 x,y=m(-3.,77.)
  47 m.plot(x,y,'b.')
  48 x,y=m(10.,77.)
  49 m.plot(x,y,'b.')
  50 title('grid points with profiles(red), spectrum (blue)')
  51 
  52 show()
  53 savefig('Map.png',dpi=75)