#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 *

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/'
crop=array([600,600,2600,4700]) # Crop region xsize,ysize,xoff,yoff
B=read_ls7(datadir,crop)
close(1)
figure(1)
display_rgb(B,crop=crop)
savefig('landsat7.png',dpi=100)

blue=B['B10']
red=B['B30']
green=B['B20']

hblue=histogram(blue,bins=256,range=[0,255],normed=True)
hred=histogram(red,bins=256,range=[0,255],normed=True)
hgreen=histogram(green,bins=256,range=[0,255],normed=True)

pdfblue, x=hblue[0],hblue[1]
pdfred, x=hred[0],hblue[1]
pdfgreen, x=hgreen[0],hblue[1]

import pickle
regions=pickle.load(open('icedata.dat','r'))

close(2)
figure(2)
plot(x,pdfblue,'b',label='blue')
plot(x,pdfred,'r',label='red')
plot(x,pdfgreen,'g',label='green')
legend()
xlabel('gray level')
ylabel('pdf')

for key in regions.keys():
	print key
	klasse=regions[key]
	blueregion=klasse['B10']	
	redregion=klasse['B30']
	greenregion=klasse['B20']
	plot([mean(blueregion)],[pdfblue[int(mean(blueregion))]],'bo')
	text(mean(blueregion),pdfblue[int(mean(blueregion))],key)
	plot([mean(redregion)],[pdfred[int(mean(redregion))]],'ro')	
	text(mean(redregion),pdfred[int(mean(redregion))],key)
	plot([mean(greenregion)],[pdfgreen[int(mean(greenregion))]],'go')
	text(mean(greenregion),pdfgreen[int(mean(greenregion))],key)
show()



}}}

[[attachment:image2.png]]