Size: 2006
Comment:
|
← Revision 5 as of 2008-12-11 10:37:10 ⇥
Size: 2074
Comment:
|
Deletions are marked like this. | Additions are marked like this. |
Line 24: | Line 24: |
B[band]=gdal.Open(fn).ReadAsArray(xsize=C[0],ysize=C[1],xoff=C[2],yoff=C[3]) | B[band]=gdal.Open(fn).ReadAsArray(ysize=C[0],xsize=C[1],yoff=C[2],xoff=C[3]) |
Line 30: | Line 30: |
C=kw['crop']#extent is (left, right, bottom, top) data values of the axes | C=kw['crop'] #extent is (left, right, bottom, top) data values of the axes |
Line 39: | Line 39: |
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 |
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 |
Line 44: | Line 44: |
Boi=read_ls7(datadir,cropoi) | Bi=read_ls7(datadir,cropi) |
Line 47: | Line 47: |
#display_rgb(B,crop=crop) | crop=array([600,600,2600,4700]) B=read_ls7(datadir,crop) display_rgb(Bw,crop=cropw) |
Line 50: | Line 55: |
Ball={'Water':Bw,'Old Ice':Boi,'Melt Ice':Bmi,'Water on Ice':Bwi} | Ball={'Water':Bw,'Ice':Bi,'Melting Ice':Bmi,'Water on Ice':Bwi} |
Line 53: | Line 58: |
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'))