GRASS Short Course - Notes for Part 4: Image Processing

Pre-processing of LANDSAT-7 data

In our example we import LANDSAT-7 thermal channel.

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",
            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 \

gdalinfo p033r029_7k20000712_z13_nn61_NAD27.tif
-> PROJCS["NAD27 / UTM zone 13N",
            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: 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 tm6_lowgain_20000712
-> min=131

# Dataset from  20000712 -> 12 June 2000, important for gain/bias selection.
# 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

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

g.region rast=tm.6 -p
r.mapcalc "tm.6rad=((17.04 - 0.)/(255. - 1.))*(tm.6 - 1.) + 0." -r tm.6rad
-> min=8.721260

# 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 temp_kelvin
-> min=294.966092

# convert to degree Celsius:
r.mapcalc "temp_celsius_12Jul2000=temp_kelvin - 273.15" -r temp_celsius_12Jul2000
-> min=21.816092

# apply new color table, display:
r.colors temp_celsius_12Jul2000 col=rules << EOF
0 blue
15 green
30 yellow
45 red
d.rast.leg temp_celsius_12Jul2000

#Remember: this are surface temperatures, not air temperatures!
#The satellite is approximately passing at 9:30 local time.

Image fusion with LANDSAT-7

Use the script

SPOT Image geocoding (rectification)

In this example we work with the Imagery location and rectify the SPOT sample data into the Spearfish location:
#Enter 'Imagery' location, then run:

#A) Rectification of SPOT PAN channel. gr=spotpan in=spot.p gr=spotpan loc=spearfish map=user1
d.mon x0

#use "roads" raster map or "spot.image" as reference (10+ GCPs):
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. gr=spotms subgr=spotms,, gr=spotms loc=spearfish map=user1

#use "roads" raster map or new "spot.prect" as reference (10+ GCPs):
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): spot.prect

Next we verify the rectified images:

#Leave GRASS, restart with 'Spearfish' location:
g.region rast=spot.prect -p
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

#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

SPOT Image classification

Now we use the rectified images to perform satellite image classification.
#Start GRASS with 'Spearfish' location:

#Unsupervised image classification: gr=spotms sub=spotms,,
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:
-> create area map 'smap.train'

#create group: gr=spotms subgr=spotms,,

#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: -r spotms.smap
d.vect.area -r spotms.smap

SPOT Image fusion

Finally we do image fusion of SPOT PAN with SPOT MS channels (download
#Brovey fusion method: pan=spot.prect out=brov

#IHS fusion method:
#use i.rgb.his/i.his.rgb
#see book for details.

See also Import of LANDSAT-7 document

Cosine correction for terrain effects to minimize shadows

The cosine method is rather poor, but we use it for illustration.
g.region n=4924413 s=4922034 w=589864 e=593032 -p
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.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 ( 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 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...

© 2003 Markus Neteler (neteler AT
Back Course HOME
Last change: $Date: 2004-02-15 20:49:53 +0100 (Sun, 15 Feb 2004) $