1 import sys,os
   2 home=os.getenv('HOME')
   3 import gdal,glob
   4 from pylab import *
   5 
   6 def sensor_params():
   7     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)}
   8     Resolution={'B10':1,'B20':1,'B30':1,'B40':1,'B50':1,'B70':1,'B61':0.5,'B62':0.5,'B80':2} # Relative to 30m
   9     return {'f':Freqs,'r':Resolution}
  10 
  11 
  12 def read_ls7(datadir,crop):
  13     SP=sensor_params()
  14     files=glob.glob(datadir+'*.TIF')
  15     B={}
  16     for fn in files:
  17         band=fn[-7:-4]
  18         C=((crop*SP['r'][band]).astype(int))
  19         print 'Reading Band ',band,'Size Position',C
  20         B[band]=gdal.Open(fn).ReadAsArray(xsize=C[0],ysize=C[1],xoff=C[2],yoff=C[3])
  21     return B
  22 
  23 def display_rgb(B,**kw):
  24     RGB=array([B['B30'],B['B20'],B['B10']]).transpose()
  25     if kw.has_key('crop'):
  26         C=kw['crop']#extent is (left, right, bottom, top) data values of the axes
  27         extent=(C[2],C[2]+C[0],C[3],C[3]+C[1])
  28         imshow(RGB,origin='lower',interpolation='nearest',extent=extent)
  29     else:
  30         imshow(RGB,origin='lower',interpolation='nearest')
  31     show()
  32     return
  33 
  34 datadir='/pf/u/u241127/satbild/LE71770012008201SGS00/'
  35 crop=array([600,600,2600,4700]) # Crop region xsize,ysize,xoff,yoff
  36 B=read_ls7(datadir,crop)
  37 close(1)
  38 figure(1)
  39 display_rgb(B,crop=crop)
  40 savefig('landsat7.png',dpi=100)
  41 
  42 blue=B['B10']
  43 red=B['B30']
  44 green=B['B20']
  45 
  46 hblue=histogram(blue,bins=256,range=[0,255],normed=True)
  47 hred=histogram(red,bins=256,range=[0,255],normed=True)
  48 hgreen=histogram(green,bins=256,range=[0,255],normed=True)
  49 
  50 pdfblue, x=hblue[0],hblue[1]
  51 pdfred, x=hred[0],hblue[1]
  52 pdfgreen, x=hgreen[0],hblue[1]
  53 
  54 import pickle
  55 regions=pickle.load(open('icedata.dat','r'))
  56 
  57 close(2)
  58 figure(2)
  59 plot(x,pdfblue,'b',label='blue')
  60 plot(x,pdfred,'r',label='red')
  61 plot(x,pdfgreen,'g',label='green')
  62 legend()
  63 xlabel('gray level')
  64 ylabel('pdf')
  65 
  66 for key in regions.keys():
  67         print key
  68         klasse=regions[key]
  69         blueregion=klasse['B10']        
  70         redregion=klasse['B30']
  71         greenregion=klasse['B20']
  72         plot([mean(blueregion)],[pdfblue[int(mean(blueregion))]],'bo')
  73         text(mean(blueregion),pdfblue[int(mean(blueregion))],key)
  74         plot([mean(redregion)],[pdfred[int(mean(redregion))]],'ro')     
  75         text(mean(redregion),pdfred[int(mean(redregion))],key)
  76         plot([mean(greenregion)],[pdfgreen[int(mean(greenregion))]],'go')
  77         text(mean(greenregion),pdfgreen[int(mean(greenregion))],key)
  78 show()