1 import sys,os
   2 home=os.getenv('HOME')
   3 import gdal,glob
   4 from pylab import *
   5 import pickle
   6 
   7 def sensor_params():
   8     Freqs={'B10':(0.45,0.52),'B20':(0.52,0.60),'B30':(0.63,0.69),'B40':(0.77,0.90),'B50':(1.55,1.75),'B70':(2.09,2.35),'B61':(10.4,12.5),'B62':(10.4,12.5),'B80':(0.52,0.90)}
   9     Resolution={'B10':1,'B20':1,'B30':1,'B40':1,'B50':1,'B70':1,'B61':0.5,'B62':0.5,'B80':2} # Relative to 30m
  10     return {'f':Freqs,'r':Resolution}
  11 
  12 
  13 def read_ls7(datadir,crop):
  14     SP=sensor_params()
  15     files=glob.glob(datadir+'*.TIF')
  16     B={}
  17     for fn in files:
  18         band=fn[-7:-4]
  19         C=((crop*SP['r'][band]).astype(int))
  20         print 'Reading Band ',band,'Size Position',C
  21         B[band]=gdal.Open(fn).ReadAsArray(ysize=C[0],xsize=C[1],yoff=C[2],xoff=C[3])
  22     return B
  23 
  24 def display_rgb(B,**kw):
  25     RGB=array([B['B30'],B['B20'],B['B10']]).transpose()
  26     if kw.has_key('crop'):
  27         C=kw['crop'] #extent is (left, right, bottom, top) data values of the axes
  28         extent=(C[2],C[2]+C[0],C[3],C[3]+C[1])
  29         imshow(RGB,origin='lower',interpolation='nearest',extent=extent)
  30     else:
  31         imshow(RGB,origin='lower',interpolation='nearest')
  32     show()
  33     return
  34 
  35 datadir='/pf/u/u241127/satbild/LE71770012008201SGS00/'
  36 cropw=array([30,35,4880,3675]) # Crop region xsize,ysize,xoff,yoff
  37 cropi=array([20,30,4735,3410]) # Crop region xsize,ysize,xoff,yoff
  38 cropmi=array([30,30,4685,3240]) # Crop region xsize,ysize,xoff,yoff
  39 cropwi=array([20,15,4865,3493]) # Crop region xsize,ysize,xoff,yoff
  40 Bw=read_ls7(datadir,cropw)
  41 Bi=read_ls7(datadir,cropi)
  42 Bmi=read_ls7(datadir,cropmi)
  43 Bwi=read_ls7(datadir,cropwi)
  44 
  45 crop=array([600,600,2600,4700])
  46 B=read_ls7(datadir,crop)
  47 
  48 
  49 display_rgb(Bw,crop=cropw)
  50 #savefig('landsat7.png',dpi=100)
  51 
  52 Ball={'Water':Bw,'Ice':Bi,'Melting Ice':Bmi,'Water on Ice':Bwi}
  53 
  54 pickle.dump(Ball,open('icedata.dat','w'))