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(xsize=C[0],ysize=C[1],xoff=C[2],yoff=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([10,20,2710,5075]) # Crop region xsize,ysize,xoff,yoff
cropoi=array([15,15,3025,5125]) # Crop region xsize,ysize,xoff,yoff
cropmi=array([10,10,3080,5140]) # Crop region xsize,ysize,xoff,yoff
cropwi=array([8,8,3097,4740]) # Crop region xsize,ysize,xoff,yoff
Bw=read_ls7(datadir,cropw)
Boi=read_ls7(datadir,cropoi)
Bmi=read_ls7(datadir,cropmi)
Bwi=read_ls7(datadir,cropwi)
#display_rgb(B,crop=crop)
#savefig('landsat7.png',dpi=100)

Ball={'Water':Bw,'Old Ice':Boi,'Melt Ice':Bmi,'Water on Ice':Bwi}

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

icedata.dat