An introduction to using Landsat data with GRASS

This brief introduction to acquiring and processing Landsat imagery is based on the notes I took while creating PAN-sharpened, false-colour mosaics of the Landsat-7 image tiles that cover my city. This was my first foray into remote sensing, image processing, and GIS, and I could have used a guide like this one when starting out.

The satellite

Landsat-7 is the latest of the Landsat satellites, which have maintained a continuous record of the Earth's surface since 1972. Launched on April 15 1999, its sun-synchronous, near-polar orbit allows it to record the surface of the Earth in a pattern of overlapping 185km swathes, completing one scan every sixteen days (or 233 orbits). The Worldwide Referencing System (WRS) catalogues this pattern as a tiled global grid of 233 paths and 248 rows, where every tile is addressed by its path and row number.

The Enhanced Thematic Mapper (ETM+) on board Landsat-7 is a multi-spectral radiometric sensor that records eight bands of data with varying spectral and spatial resolutions (30m spatial resolution for red, green, blue, near infrared, and two bands of medium infrared; 60m for thermal infrared; and a 15m panchromatic band).

Each raw scene obtained from Landsat-7 undergoes levels of processing to remove radiometric and geometric errors. (The details of this process are beyond the scope of this introduction. I shall assume that you are using corrected data from one of the sources listed below.)

The data

The Global Land Cover Facility provides free Web/FTP access to the orthorectified GeoCover data from Landsat. This data contains at least one scene (and often two or more) from every part of the Earth's surface except Antarctica, and it may be used without any restrictions for research, educational, or commercial purposes. (But please take note of the request for proper citation). The GLCF also provides access to 90m SRTM DEM and other useful data. My appreciation for them cannot be overstated.

NASA has a list of other sources (none of which I have used personally).

The GeoCover data for each tile is usually available as a set of eight GeoTIFF files (one for each ETM+ band). The bands, their characteristics, and examples of the corresponding file names are as shown below.

Thus, a file named "p146r040_7t19991022_z43_nn10.tif" is the scene from path 146/row 040, collected by Landsat-7's t sensor (30m multi-spectral) on 1999-10-22, projected into the UTM zone 43, for band 1. (Note that the data for the thermal band #6 is collected at both low (nn61) and high (nn62) gain settings.)

Data from the various bands can be combined to form true-colour (RGB) or false-colour (G, NIR, MIR) composites. A composite of some 30m bands can be resampled and combined with the 15m panchromatic band data to yield a high-resolution, coloured image. Two or more adjacent overlapping scenes may be combined into a mosaic, and so on. I wanted to create a sharpened mosaic of two adjacent false-colour composite tiles that covered Delhi.

GRASS

I used the open-source GRASS GIS to create the composites and mosaic.

A useful introduction to GRASS would deserve its own web page, which I shall not attempt to write. If you have not used GRASS before, I must, for now, refer you to the documents that helped me to get started: A quick start guide for GRASS 6, an old article introducing GRASS, and an assortment of other tutorials.

I used r.in.gdal to import each TIFF band into GRASS as a raster map, creating a new location in the process.

# Create a new location with the first map
r.in.gdal input=p146r040_7t19991022_z43_nn20.tif output=ls_146_40_20 \
          location=landsat

# Import the other bands into that location
g.mapset location=landsat mapset=PERMANENT
r.in.gdal input=p146r040_7t19991022_z43_nn40.tif output=ls_146_40_40
r.in.gdal input=p146r040_7t19991022_z43_nn50.tif output=ls_146_40_50
r.in.gdal input=p146r040_7p19991022_z43_nn80.tif output=ls_146_40_80

I used i.fusion.brovey to create a composite of bands 2, 4, and 5 (green, near infrared, and medium infrared), sharpened with the data from band 8 (panchromatic); and combined the three RGB output maps into one with r.composite. The Brovey transformation apparently does better with the infrared channels than while producing true-colour composites.

# Create ls2_146_40.{red,green,blue} maps
i.fusion.brovey -l ms1=ls_146_40_20 ms2=ls_146_40_40 \
                   ms3=ls_146_40_50 pan=ls_146_40_80 \
                   outputprefix=ls2_146_40

# Notice the reversal of the R/G channels
r.composite g_map=ls2_146_40.red r_map=ls2_146_40.green \
            b_map=ls2_146_40.blue output=ls2_146_40_rgb

I also created a PAN-sharpened ls2_147_40 raster map as shown above, and used r.null (to remove the black image borders) and r.patch to combine the two scenes into an overlapping mosaic.

g.region rast=ls2_146_40_rgb,ls2_147_40_rgb

r.null map=ls2_147_40_rgb setnull=0
r.null map=ls2_146_40_rgb setnull=0

r.patch in=ls2_147_40_rgb,ls2_146_40_rgb out=ls2_14x_40_rgb

The results — especially with IR composites — are excellent. The output map may be viewed with d.rast, or exported into a file with with an r.out.* command like r.out.png.

d.mon start=x0
d.rast ls2_14x_40_rgb

r.out.png input=ls2_14x_40_rgb output=mosaic.png

Alternatives to GRASS

Unfortunately, GRASS is not as easy to set up or work with as could be hoped for. Furthermore, i.fusion.brovey is slow (sometimes taking more than ninety minutes to create a 17236*15216 PAN-sharpened composite for a single scene, on a 3GHz Pentium 4 with 2GB of RAM.) I have not found a satisfactory alternative yet, but this section discusses two of the alternatives I've tried: ImageMagick and netpbm.

ImageMagick

composite -compose CopyGreen p147r040_7t20000913_z43_nn50.tif \
          p147r040_7t20000913_z43_nn40.tif rg.pnm
composite -compose CopyBlue p147r040_7t20000913_z43_nn20.tif \
          rg.pnm rgb.pnm
convert -sample 17094x15224 rgb.pnm rgb2.pnm
composite rgb2.pnm p147r040_7p20000913_z43_nn80.tif rgbp.pnm

This is much faster than GRASS (especially if the intermediate files are PNM), finishing in less than 20 minutes on my machine. Sadly, the result is not nearly as good. The colours are muddy, and the image less sharp overall.

netpbm

tifftopnm p147r040_7t20000913_z43_nn20.tif > 20.pgm
tifftopnm p147r040_7t20000913_z43_nn50.tif > 50.pgm
tifftopnm p147r040_7t20000913_z43_nn40.tif > 40.pgm
tifftopnm p147r040_7p20000913_z43_nn80.tif > 80.pgm

rgb3toppm 40.pgm 50.pgm 20.pgm > rgb.ppm
pnmscale 2 rgb.ppm > rgb2.ppm
pnmnorm -bvalue 18 -wvalue 156 > pan.pgm < 80.pgm
pnmarith -add pan.pgm rgb2.ppm | pnmgamma 0.66 > rgbp.ppm

This is both faster (~10 minutes) and better than using ImageMagick, and the output is clear and sharp, but it's still not quite as good as those obtained with GRASS, and the performance is not always good (using it to composite the visible bands is distinctly less useful). Still, it may be possible to tweak the parameters to obtain a better result. I haven't tried too hard.

GeoTIFF files are sometimes encoded with 32 bits per sample, and many programs can't read such files. One way to convert them into something usable is "gdal_translate -ot Uint16 -of PNM i.tiff o.pnm". A related problem is that most image processing programs can't copy the GeoTIFF referencing tags into the output, so that data is lost. No doubt GDAL could be used to copy this information (as this python script does), but I haven't investigated further.

There are commercial programs like AlphaPixel's PixelSense that composite Landsat images. I haven't tried any of them (they're not available on Linux, and I couldn't afford them if they were).

Suggestions with respect to this section are very welcome.

References

• FAQs: Landsat-7, Satellite Imagery (old), GLCF, GRASS.
• Remote Sensing Tutorial (and more tutorials).
• FreeGIS directory.
• GRASS tutorials for 5.7, an older version, and another that is older still, and poorly formatted to boot.

Please send questions, corrections, and suggestions to ams@toroid.org.