#acl AdminGroup:read,write,delete,revert  EditorGroup:read,write StudentGroup:read,write,delete,revert All:read
#format wiki
{{{#!python
import sys,os
home=os.getenv('HOME')
import gdal,glob
from pylab import *
import pickle

def sensor_params():
    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)}
    Resolution={'B10':1,'B20':1,'B30':1,'B40':1,'B50':1,'B70':1,'B61':0.5,'B62':0.5,'B80':2} # Relative to 30m
    return {'f':Freqs,'r':Resolution}


def read_ls7(datadir,crop):
    SP=sensor_params()
    files=glob.glob(datadir+'*.TIF')
    B={}
    for fn in files:
        band=fn[-7:-4]
        C=((crop*SP['r'][band]).astype(int))
        print 'Reading Band ',band,'Size Position',C
        B[band]=gdal.Open(fn).ReadAsArray(ysize=C[0],xsize=C[1],yoff=C[2],xoff=C[3])
    return B

def display_rgb(B,**kw):
    RGB=array([B['B30'],B['B20'],B['B10']]).transpose()
    if kw.has_key('crop'):
        C=kw['crop'] #extent is (left, right, bottom, top) data values of the axes
        extent=(C[2],C[2]+C[0],C[3],C[3]+C[1])
        imshow(RGB,origin='lower',interpolation='nearest',extent=extent)
    else:
        imshow(RGB,origin='lower',interpolation='nearest')
    show()
    return

datadir='/pf/u/u241127/satbild/LE71770012008201SGS00/'
cropw=array([30,35,4880,3675]) # Crop region xsize,ysize,xoff,yoff
cropi=array([20,30,4735,3410]) # Crop region xsize,ysize,xoff,yoff
cropmi=array([30,30,4685,3240]) # Crop region xsize,ysize,xoff,yoff
cropwi=array([20,15,4865,3493]) # Crop region xsize,ysize,xoff,yoff
Bw=read_ls7(datadir,cropw)
Bi=read_ls7(datadir,cropi)
Bmi=read_ls7(datadir,cropmi)
Bwi=read_ls7(datadir,cropwi)

crop=array([600,600,2600,4700])
B=read_ls7(datadir,crop)


display_rgb(Bw,crop=cropw)
#savefig('landsat7.png',dpi=100)

Ball={'Water':Bw,'Ice':Bi,'Melting Ice':Bmi,'Water on Ice':Bwi}

pickle.dump(Ball,open('icedata.dat','w'))


}}}

[[attachment:icedata.dat]]