In this second post of the Landsat 8 series, we have a look at how to create simple color composites with uncorrected data, a common and useful task that is often misunderstood by beginners. For this purpose, it is necessary to understand two concepts: how digital sensors register information and how color images are constructed.

The first key concept is that sensors on board satellites or air planes perform radiometric measurements in specific ranges of wavelength of the electromagnetic spectrum. Sensor elements capture a portion of the outgoing (from Earth) radiation in a given spectral window, which is then converted to digital numbers, stored and, together with the set of neighboring measurements, coded as regular image grids. Images acquired in specific windows of the spectrum are called ‘bands’. The spectral ranges in which they were acquired, measured in nanometers or micrometers, together with their nominal pixel size, are key properties of the sensors. Then, each band will be an image file. For Landsat 8 we have:

Designation Range (µm) Pixel size Denomination/Use
Band 1 0.43 - 0.45 30 Coastal aerosol
Band 2 0.45 - 0.51 30 Blue
Band 3 0.53 - 0.59 30 Green
Band 4 0.64 - 0.67 30 Red
Band 5 0.85 - 0.88 30 Near Infrared (NIR)
Band 6 1.57 - 1.65 30 Short wave infrared (SWIR) 1
Band 7 2.11 - 2.29 30 Short wave infrared (SWIR) 2
Band 8 0.50 - 0.68 15 Panchromatic
Band 9 1.36 - 1.38 30 Cirrus
Band 10 10.60 - 11.19 100 * (30) Thermal Infrared (TIRS) 1
Band 11 11.50 - 12.51 100 * (30) Thermal Infrared (TIRS) 2

The second key concept is that color images can be formed by projecting three separate gray-level (black and white) images through three separate color channels: blue, green, red (in increasing order of wavelength). This is what television sets and projectors do. If the three images are exactly the same, the result will be another gray-level image, but if these images contain differences, these will become manifest as colored areas, depending on the relative magnitude of the levels at a given point.

In the specific case of satellite images, we can create a ‘true-color’ image by combining the three images acquired in the blue, green and red parts of the visible part of the spectrum, through the blue, green and red channels of the image. For example, to create a true color composite of the Bay of Naples, we put the image collected in the red part of the visible spectrum (band 4 in Landsat 8) in the red channel, the green image (band 3) on the green channel, and the blue image (band 2) on the blue channel. We can do this from the command line by simply using the GDAL tool gdal_merge.py:

> gdal_merge.py -separate -o naples_true_composite.tif LC81890322013210LGN01_B4.TIF LC81890322013210LGN01_B3.TIF LC81890322013210LGN01_B2.TIF
0...10...20...30...40...50...60...70...80...90...100 - done.

Our test image will look like this:

Bay of Naples, Italy
Bay of Naples, Italy

Landsat images are about 185×185 km in size. We might not need such a big area, and we then might need to clip a subset. To achieve this from the command line, for example, to extract the Vesuvius area, we can use the GDAL tool gdal_translate:

> gdal_translate -projwin 431000 4539000 471000 4489000 naples_true_composite.tif vesuvius.tif
Input file size is 7761, 7901
Computed -srcwin 2317 1460 1333 1667 from projected window.
0...10...20...30...40...50...60...70...80...90...100 - done.

where, using the option -projwin we introduced the window in projected coordinates: ulx (upper left x), uly (upper left y), lrx (lower right x), lry (lower right y). Note that we will need to pick the coordinates of the corners. We can do this in any viewer, such as QGIS.

Our image now looks like this:

Vesuvius area, Campania, Italy.
Vesuvius area, Campania, Italy.

which is basically an approximately 20×20 km area around Vesuvius. Note the northern tip of Capri on the southwest corner.

Alternatively, we might create a ‘false color composite’, by using a different combination of bands in different channels. An extremely useful combination, particularly for vegetation studies, is to use the infrared band, which is not visible to the human eye, in the red channel, the red image in the green channel, and the green image in the blue channel (remember that images are usually specified in order: RGB). Using the same procedures as above, we do:

> gdal_merge.py -separate -o naples_false_composite.tif LC81890322013210LGN01_B5.TIF LC81890322013210LGN01_B4.TIF LC81890322013210LGN01_B3.TIF
0...10...20...30...40...50...60...70...80...90...100 - done.

and we get:

Bay of Naples, Italy, infrared false color composite.
Bay of Naples, Italy, infrared false color composite.

Note the strong red color of the vegetation, as it is much more reflective in the infrared than in the visible wavelengths.

As a final point. We can check the properties of the Vesuvius image with gdalinfo:

> gdalinfo vesuvius.tif
Driver: GTiff/GeoTIFF
Files: vesuvius.tif
Size is 1333, 1667
Coordinate System is:
PROJCS["WGS 84 / UTM zone 33N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",15],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32633"]]
Origin = (430995.000000000000000,4539015.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 430995.000, 4539015.000) ( 14d10'46.25"E, 40d59'57.85"N)
Lower Left ( 430995.000, 4489005.000) ( 14d11' 6.14"E, 40d32'56.15"N)
Upper Right ( 470985.000, 4539015.000) ( 14d39'17.97"E, 41d 0' 6.51"N)
Lower Right ( 470985.000, 4489005.000) ( 14d39'26.33"E, 40d33' 4.67"N)
Center ( 450990.000, 4514010.000) ( 14d25' 9.19"E, 40d46'32.18"N)
Band 1 Block=1333x1 Type=UInt16, ColorInterp=Gray
Band 2 Block=1333x1 Type=UInt16, ColorInterp=Undefined
Band 3 Block=1333x1 Type=UInt16, ColorInterp=Undefined

We see that the image has a completely defined spatial reference system (SRS), with datum, coordinates referred to the WGS84 ellipsoid, and a cartographic projection: Universal Transverse Mercator (UTM) zone 33 north. We also see from the last lines that this image is coded as in 16-bit unsigned integers. In the next post, we will see how to convert these digital numbers into geophysical measurements of radiance.

Advertisements