There are two gain states (low and high gain). The science goal in switching gain states is to maximize the instrument's 8 bit radiometric resolution without saturating the detectors. For thermal data, both low and high gain data are available by default. For other bands (1 to 5, and 7) the satellite will acquire image data in one of two possible gain settings. The low gain setting measures a greater radiance but with decreased sensor sensitivity (for very bright regions to avoid detector saturation), while high gain measures a lesser radiance range but with increased sensitivity (over areas of low reflectance). See also here.
First we reproject the data to our target projection/datum/ellipsoid, then we cut out the area of interest:
gdalinfo p033r029_7k20000712_z13_nn61.tif
-> PROJCS["WGS 84 / UTM zone 13N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.2572235629972,
[...]
#Since Spearfish Location is in UTM/NAD27/Clarke66, we have to reproject:
#EPSG file: /usr/local/share/proj/epsg
gdalwarp -t_srs '+init=epsg:26713' p033r029_7k20000712_z13_nn61.tif \
p033r029_7k20000712_z13_nn61_NAD27.tif
#validate:
gdalinfo p033r029_7k20000712_z13_nn61_NAD27.tif
-> PROJCS["NAD27 / UTM zone 13N",
GEOGCS["NAD27",
DATUM["North_American_Datum_1927",
SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
[...]
#now ready to import into GRASS/Spearfish Location. But we only need
#a small portion of that area. So we cut it out:
grass5 ~/grassdata/spearfish/user1/
g.region -dp
projection: 1 (UTM)
zone: 13
datum: nad27
ellipsoid: clark66
north: 4928010
south: 4913700
west: 589980
east: 609000
nsres: 30
ewres: 30
rows: 477
cols: 634
#boundary coordinates: W N E S
gdal_translate -of GTiff -projwin 589980 4928010 609000 4913700 \
p033r029_7k20000712_z13_nn61_NAD27.tif p033r029_7k20000712_z13_nn61_NAD27_small.tif
#import:
r.in.gdal in=p033r029_7k20000712_z13_nn61_NAD27_small.tif out=tm6_lowgain_20000712
d.rast tm6_lowgain_20000712
d.vect roads
# (if 'roads' vector map is not there, run 'v.convert' first)
The LANDSAT-7 low-gain temperature channel is imported now. Next step is the conversion of the pixel values to real data (since they are coded). Here we want to convert to Kelvin and later to degree Celsius.
r.info -r tm6_lowgain_20000712 -> min=131 max=175 # Dataset from 20000712 -> 12 June 2000, important for gain/bias selection. # http://geog.tamu.edu/~liu/courses/g361/note9.pdf # http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_htmls/chapter11/chapter11.html # Calculate DN to spectral radiances: # radiance = gain * DN + offset # When using GeoTIFF 1G products, only LMIN, LMAX, QCALMIN and QCALMAX are given. # Convert TM7/b61 digital numbers (DN) to spectral radiances (apparent radiance at sensor) # radiance = gain * DN + offset # = ((LMAX-LMIX)/(QCALMAX-QCALMIN))*(DN-QCALMIN)+LMIN grep BAND61 p033r029_7x20000712.met # BAND61_FILE_NAME = "p033r029_7k20000712_z13_nn61.tif" # LMAX_BAND61 = 17.040 # LMIN_BAND61 = 0.000 # QCALMAX_BAND61 = 255.0 # QCALMIN_BAND61 = 1.0 # CORRECTION_METHOD_GAIN_BAND61 = "CPF" # STRIPING_BAND61 = "NONE" g.region rast=tm.6 -p r.mapcalc "tm.6rad=((17.04 - 0.)/(255. - 1.))*(tm.6 - 1.) + 0." r.info -r tm.6rad -> min=8.721260 max=11.673071 # convert spectral radiances to absolute temperatures: # T = K2/ln(K1/L_l + 1)) r.mapcalc "temp_kelvin=1260.56/(log (607.76/tm.6rad + 1))" r.info -r temp_kelvin -> min=294.966092 max=315.820713 # convert to degree Celsius: r.mapcalc "temp_celsius_12Jul2000=temp_kelvin - 273.15" r.info -r temp_celsius_12Jul2000 -> min=21.816092 max=42.670713 # apply new color table, display: r.colors temp_celsius_12Jul2000 col=rules << EOF 0 blue 15 green 30 yellow 45 red EOF d.rast.leg temp_celsius_12Jul2000 #Remember: this are surface temperatures, not air temperatures! #The satellite is approximately passing at 9:30 local time.
#Enter 'Imagery' location, then run: #A) Rectification of SPOT PAN channel. i.group gr=spotpan in=spot.p i.target gr=spotpan loc=spearfish map=user1 d.mon x0 #use "roads" raster map or "spot.image" as reference (10+ GCPs): i.points i.rectify gr=spotpan ext=rect order=3 #... after a while you receive an email that the image(s) are rectified. #B) Rectification of SPOT Green + Red + NIR (MS) channels. i.group gr=spotms subgr=spotms in=spot.ms.1,spot.ms.2,spot.ms.3 i.target gr=spotms loc=spearfish map=user1 #use "roads" raster map or new "spot.prect" as reference (10+ GCPs): i.points i.rectify gr=spotms ext=rect order=3 #... after a while you receive an email that the image(s) are rectified. #Note down the acquisition date of the SPOT PAN image for later (terrain #effects correction): r.info spot.prect
Next we verify the rectified images:
#Leave GRASS, restart with 'Spearfish' location: g.region rast=spot.prect -p d.erase d.rast spot.prect #apply better color table: r.colors spot.prect col=grey.eq d.rast spot.prect #verify with other map: d.vect roads col=red d.zoom #check if the accuracy is sufficient. otherwise go back to the 'Imagery' #location and add/replace GCPs in i.points. It remembers the GCPs already #defined.
#Start GRASS with 'Spearfish' location: #Unsupervised image classification: i.group gr=spotms sub=spotms in=spot.ms.1rect,spot.ms.2rect,spot.ms.3rect i.cluster gr=spotms sub=spotms sig=cluster class=10 i.maxlik gr=spotms sub=spotms sig=cluster class=spot.class rej=spot.rej d.rast.leg spot.class #Supervised image classification with SMAP: #digitize training areas: r.digit -> create area map 'smap.train' #create group: i.group gr=spotms subgr=spotms in=spot.ms.1,spot.ms.2,spot.ms.3 #calculate statistics: i.gensigset train=smap.train gr=spotms subgr=spotms sig=smap #apply statistics to map (radiometric/geometric SMAP approach): i.smap gr=spotms subgr=spotms sig=smap out=spotms.smap d.rast spotms.smap #vectorization of resulting polygons: r.poly -l in=spotms.smap out=spotms.smap #in 5.3: v.support -r spotms.smap d.vect.area -r spotms.smap
#Brovey fusion method: i.brovey.spot spot1=spot.ms.1rect spot2=spot.ms.2rect spot3=spot.ms.3rect pan=spot.prect out=brov #IHS fusion method: #use i.rgb.his/i.his.rgb #see book for details.
g.region n=4924413 s=4922034 w=589864 e=593032 -p
d.erase
d.rast spot.prect
#interpolation of high-res DEM (to 3m resolution):
g.region res=3 -p
#bilinear interpolation (NOTE: better to use RST interpolation!):
r.bilinear elevation.dem out=elev.bilin
#due to bug in r.bilinear we must set 0 -> NULL
r.null elev.bilin setnull=0
#copy over original color table:
r.colors elev.bilin rast=elevation.dem
#zoom to new map:
g.region rast=elev.bilin -p
d.erase
d.rast elev.bilin
d.rast spot.prect
#now calculate slope/aspect
r.slope.aspect elev.bilin as=as3 sl=sl3
#shadow correction:
#get sun position for SPOT PAN image acquisition date (r.info in Imagery
#location, see above):
r.sunmask -s y=1989 mon=5 d=27 h=10 min=58 s=50 timez=-7 elev=elevation.dem out=dummy
Using map center coordinates
Calculating sun position... (using solpos (V. 11 April 2001) from NREL)
1989.05.27, daynum 147, time: 10:58:50
long: -103.851006, lat: 44.458466, timezone: -7.000000
Solar position: sun azimuth 149.979584,
sun angle above horz.(refraction corrected) 64.404373
Sunrise time (without refraction): 04:22
Sunset time (without refraction): 19:23
No map calculation requested. Finished.
#-> we need sun angle above horz and sun azimuth
#generate new aspect map oriented from North, orientation clockwise:
#we are using an internal variable 'as':
r.mapcalc "as3_north=eval(as=360. - as3 + 90., if(as > 360., as - 360., as))"
d.rast as3_north
#generate incidence map from sun position:
r.mapcalc "cos_i=cos(90.-64.4)*cos(sl3) + sin(90.-64.4) * sin(sl3) * cos(150. - as3_north)"
r.info -r cos_i
r.colors cos_i col=grey
d.rast cos_i
d.vect roads col=red
#compare cos_i map with SPOT PAN shadows:
d.rast spot.prect
#generate contour lines from elev10:
r.contour -q elev.bilin step=20 out=contour10
d.vect contour10 col=green
d.vect.labels -m contour10 col=green att=cat
#apply the cosine correction:
r.mapcalc "factor=float(cos(90.-64.4)/cos_i)"
r.mapcalc "factor_filt=if(factor > -5.0 && factor < 5.0, factor, null())"
r.mapcalc "spot.pcorr=spot.prect * factor_filt"
r.colors spot.pcorr col=aspect
d.frame -e
d.frame -c at="0,100,50,100"
d.frame -c at="0,100,0,50"
d.rast spot.prect
echo "Original SPOT-1 PAN" | d.text col=white
#select right frame by clicking:
d.frame -s
d.rast spot.pcorr
echo "Terrain corrected SPOT-1 PAN" | d.text col=white
#-> looks acceptable with the limits of the cosine correction (overshoots)
#to erase the screen with frames:
d.frame -e
###########
#NOTE better to use 'Minnaert correction' or other methods...