The corrected image spot.pcorr has reduced topographic effects as shown in Figure 9.9 for a subregion south of Spearfish. Note that cast shadows resulting from relief (mountainous regions) are not taken into account. This approach is only a local correction. It may be extended with cast shadow maps generated by r.sunmask. Another simple method to reduce terrain effects is channel rationing which is discussed below.

9.4.3 Application: Deriving a surface temperature map from thermal channel

Several satellites such as ASTER/TERRA, LANDSAT-TM5 and LANDSAT-TM7 provide thermal channels. The data delivered by a thermal channel (channel 6 for LANDSAT systems) can be calibrated to a surface temperature map. These surface temperatures must not be confused with air temperatures. Note that the methods are different for LANDSAT-TM5 and LANDSAT-TM7, as their gain/bias values are different. For an absolute calibration of satellite-derived temperatures, atmospheric correction has to be taken into account.

Surface temperature map from LANDSAT-TM5 channel 6. The following calculations derive the effective at-satellite temperatures (LANDSAT-TM5) of the viewed earth-atmosphere system under an assumption of unity emissivity and using pre-launch calibration constants. In a first step, the gain/bias values are applied to the thermal channel to receive spectral radiances. Then the resulting pixel values are converted to absolute temperature in Kelvin. Optionally it can be recalculated to a degree Celsius temperature map:

r.mapcalc factor=float(cos(90.-64.4)/cos i) #filter out outlier:

r.mapcalc factor filt=if(factor>-5.0 && factor<5.0,factor,1)

#apply the cosine correction:

r.mapcalc spot.pcorr=spot.p * 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.p

echo Original SPOT-1 PAN I d.text col=white

#select right frame by clicking: d.frame -s d.rast spot.pcorr

echo Terrain corrected SPOT-1 PAN I d.text col=white

The map tm.temp celsius shows the distributed emitted thermal radiation in degree Celsius. The surface brightness temperature is the actual surface temperature only when the emissivity of the object in a particular waveband equals to 1.0. For most surfaces, where the emissivity is near, but not equal to 1.0, a calibration according to the Stefan-Boltzmann equation is needed when interpreting the results. You may generate a land use/land cover map through image reclassification as shown later in Section 9.8 where you apply the individual emissivity factors according to land use. With r.mapcalc you can calibrate the landuse corrected temperature map from these maps. A modified formula for deriving LANDSAT-TM5 surface temperatures was proposed by Singh, 1988 and other authors.

Surface temperature map from LANDSAT-TM7 channel 6. As in the previous case, the LANDSAT-TM7 image data have to be converted from digital numbers to spectral radiances by applying the gain/bias values. Depending on the data format, these values may be retrieved from the image metadata with gdalinfo In our example we use the LANDSAT-TM7 scene as prepared in Section 3.3.3 for the Spearfish sample LOCATION. However, the data are provided in GeoTIFF format. This requires to look up the gain and bias parameters from the accompanying metadata file which is separared from the image data.

We use the low gain thermal channel tm7 2 000071.6 as imported earlier (see Section 9.2.1). The conversion procedure is:7

# convert TM5/b6 digital numbers (DN) to spectral radiances

# (apparent radiance at sensor): radiance = gain * DN + offset g.region rast=tm.6 -p

r.mapcalc tm.6rad=l.238+(15.600-1.238)* tm.6 /255.

# convert spectral radiances to absolute temperatures:

# T = K2/ln(Kl/L l + 1))

r.mapcalc tm.temp kelvin=1260.56/(log (607.76/tm.6rad + 1.))

# convert to degree Celsius:

r.mapcalc tm.temp celsius=tm.temp kelvin - 273.15

# apply new color table, display: r.colors tm.temp celsius col=rules << EOF 0 blue

15 green 30 yellow 45 red EOF

d.rast tm.temp celsius d.legend tm.temp celsius

#rename for convenience g.rename rast=tm7 2000071.6,tra.6

grep BAND61 p033r029 7x20000712.met

# BAND61 FILE NAME = p033r029 7k20000712 zl3 nn61.tif

# LMAX BAND61 = 17.040

# LMIN BAND61 = 0.000

# QCALMAX BAND61 = 255.0

# QCALMIN BAND61 =1.0

# C0RRECTI0N METH0D GAINJAND61 = CPF

# STRIPING BAND61 = NONE

# convert TM7/b61 digital numbers (DN) to spectral radiances

# (apparent radiance at sensor)

# radiance = gain * DN + offset

# = ((LMAX-LMIX)/(QCALMAX-QCALMIN))*(DN-QCALMIN)+LMIN g.region rast=tm.6 -p

r.mapcalc tm.6rad=( (17.04 - 0.)/(2S5. - l.))*(tm.6 - 1.) + 0.

#convert spectral radiances to absolute temperatures:

# T = К2/1П(Kl/L l + 1))

# Kl: 666.09 W/ln(m 2 * sr * urn)

# K2: 1282.71 Kelvin

r.mapcalc tm.temp kelvin=1282.71 /log(666.09 / tm.6rad + 1.) #calculate degree Celsius:

r.mapcalc tm.temp celsius=tm.temp kelvin - 273.15

#apply new color table, display;

r.colors tm.temp celsius col=rules << EOF

0 blue

15 green

30 yellow

45 red

#display the map, overlay to ETM/PAN (B80) 15m channel:

g.region rast=etmpan -p

d.his i=etmpan h=tm.temp celsius

d.legend tm.temp celsius

The resulting temperature map (in degree Celsius) represents the uncorrected surface temperatures at image acquisition time (around 9:30h), see notes above for emissivity correction. For deriving these maps from other satellites such as ASTER/TERRA, please refer to the related documents.8