|
Size: 2012
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 16: | Line 16: |
| def read_ls7(datadir): | def read_ls7(datadir,crop): |
| Line 22: | Line 22: |
| #C=((crop*SP['r'][band]).astype(int)) print 'Reading Band ',band, B[band]=gdal.Open(fn).ReadAsArray() |
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]) |
| 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([100,150,2600,4900]) # Crop region ysize,xsize,yoff,xoff #cropoi=array([30,20,3085,5155]) # Crop region xsize,ysize,xoff,yoff #cropmi=array([15,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) |
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) |
| Line 48: | Line 48: |
| #crop=array([600,600,2600,4700]) B=read_ls7(datadir) figure(9) display_rgb(B) savefig('landsat7.png',dpi=100) |
crop=array([600,600,2600,4700]) B=read_ls7(datadir,crop) |
| Line 54: | Line 51: |
| #Ball={'Water':Bw,'Old Ice':Boi,'Melt Ice':Bmi,'Water on Ice':Bwi} | |
| Line 56: | Line 52: |
| #pickle.dump(Ball,open('icedata.dat','w')) | 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')) |
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'))