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 
  12 filename='AERONET_Hornsund.2009072.aqua.250m.jpg'
  13 filename='AERONET_Hornsund.2009072.terra.250m.jpg'
  14 
  15 a=load_img(filename)
  16 
  17 a=nd.zoom(a[:,:,0],0.1)
  18 
  19 north=80.2363
  20 south=73.7612
  21 east=40.4156
  22 west=-9.3026
  23 
  24 filename='AERONET_Hornsund.2009072.terra.250m.jgw'
  25 cds=array(open(filename).read().split('\n')[:-1]).astype(float)
  26 
  27 #north=cds[5]
  28 #west=cds[4]
  29 
  30 
  31 [yn,xn]=shape(a)
  32 
  33 # make grid
  34 x=linspace(west,east,xn)
  35 y=linspace(north,south,yn)
  36 
  37 [lons,lats] = meshgrid(x,y)
  38 
  39 
  40 
  41 figure(2)
  42 m = Basemap(width=2400000,height=1600000,projection='stere',lat_ts=77.0,lon_0=15.0,lat_0=77.0,resolution='l')
  43 
  44 
  45 xm,ym=m(lons,lats)
  46 
  47 xi=linspace(m.llcrnrx,m.urcrnrx,xn) # define the new grid
  48 yi=linspace(m.llcrnry,m.urcrnry,yn)
  49 modis_img_nmpg=griddata(xm.flatten(),ym.flatten(),a.flatten(),xi,yi)
  50 
  51 
  52 #m.imshow(xm,ym,a[:,:,0]) # geht nicht!
  53 m.imshow(modis_img_nmpg,cm.bone, interpolation='bilinear',aspect='auto')
  54 m.drawcoastlines()
  55 m.fillcontinents(color='gray',lake_color='aqua')
  56 m.drawmapboundary(fill_color='aqua')
  57 
  58 m.drawmeridians(np.arange(0,360,5),labels=[1,0,0,1])
  59 m.drawparallels(np.arange(-90,90,5),labels=[0,1,0,1])
  60 
  61 show()
  62 savefig('modis_test.png')