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