1 import Image
   2 from pylab import *
   3 from mpl_toolkits.basemap import Basemap
   4 import numpy as np
   5 
   6 def load_img(filename):
   7     im=Image.open(filename)
   8     return resize(fromstring(im.tostring(),uint8),(im.size[1],im.size[0],3))
   9 
  10 import scipy.ndimage as nd
  11 # define coordinates of profile
  12 lat_glatt=[77.0,78.0,79.0]
  13 lon_glatt_start,lon_glatt_stop=2.0,10.0
  14 
  15 
  16 
  17 filename='AERONET_Hornsund.2009072.aqua.250m.jpg'
  18 
  19 
  20 a=load_img(filename)
  21 
  22 a=nd.zoom(a[:,:,0],0.2)
  23 
  24 north=80.2363
  25 south=73.7612
  26 east=40.4156
  27 west=-9.3026
  28 
  29 filename='AERONET_Hornsund.2009072.aqua.250m.jgw'
  30 cds=array(open(filename).read().split('\n')[:-1]).astype(float)
  31 
  32 
  33 
  34 [yn,xn]=shape(a)
  35 
  36 # make grid
  37 x=linspace(west,east,xn)
  38 y=linspace(north,south,yn)
  39 
  40 [lons,lats] = meshgrid(x,y)
  41 
  42 start=find(x>lon_glatt_start)[0]
  43 stop=find(x>lon_glatt_stop)[0]
  44 
  45 line=[]
  46 for i in lat_glatt:
  47     line.append(find(y<i)[0])
  48 
  49 figure(2)
  50 
  51 m=Basemap(projection='merc',llcrnrlat=74,urcrnrlat=80,llcrnrlon=-9,urcrnrlon=20,resolution='h')
  52 
  53 xm,ym=m(lons,lats)
  54 
  55 xi=linspace(m.llcrnrx,m.urcrnrx,xn) # define the new grid
  56 yi=linspace(m.llcrnry,m.urcrnry,yn)
  57 modis_img_nmpg=griddata(xm.flatten(),ym.flatten(),a.flatten(),xi,yi)
  58 
  59 
  60 m.imshow(modis_img_nmpg,cm.bone, interpolation='bilinear',aspect='auto')
  61 m.drawcoastlines()
  62 m.fillcontinents(color='gray',lake_color='aqua')
  63 m.drawmapboundary(fill_color='aqua')
  64 
  65 m.drawmeridians(np.arange(0,360,5),labels=[1,0,0,1])
  66 m.drawparallels(np.arange(-90,90,2),labels=[0,1,0,1])
  67 
  68 cls=('ryg')
  69 for k,i in enumerate(line):
  70     x1,y1=m(x[start],y[i])
  71     x2,y2=m(x[stop],y[i])
  72     m.plot([x1,x2],[y1,y2],'-'+cls[k],linewidth=3)
  73 
  74 
  75 
  76 for j in range(182,191):
  77     x,y=m(j-180.,77)
  78     m.plot(x,y,'b.',markersize=12)
  79     x,y=m(j-180.,78)
  80     m.plot(x,y,'b.',markersize=12)
  81     x,y=m(j-180.,79)
  82     m.plot(x,y,'b.',markersize=12)    
  83 
  84 show()
  85 savefig('modis_test.png')

   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 plot(Temp[:,13,187],Lv)
  24 ylim((1000.0,1.0))
  25 xlim((200.0,280.0))
  26 xlabel('T in Kelvin')
  27 ylabel('p in hPa')
  28 title('used profiles lon=77')
  29 show()
  30 savefig('Soundings1.png',dpi=75)
  31 figure(2)
  32 plot(Temp[:,12,184],Lv)
  33 ylim((1000.0,1.0))
  34 xlim((200.0,280.))
  35 xlabel('T in Kelvin')
  36 ylabel('p in hPa')
  37 title('used profiles lon=78')
  38 show()
  39 savefig('Soundings2.png',dpi=75)
  40 figure(3)
  41 plot(Temp[:,11,187],Lv)
  42 ylim((1000.0,1.0))
  43 xlim((200.0,280.0))
  44 xlabel('T in Kelvin')
  45 ylabel('p in hPa')
  46 title('used profiles lon=79')
  47 show()
  48 savefig('Soundings3.png',dpi=75)