In the last post of this series we will review the necessary steps for performing topographic and illumination correction of the MODIS images. Topographic and illumination corrections are important when we need to obtain the apparent surface reflectance as an end product. This is necessary when we want to identify materials on the Earth’s surface by deriving empirical spectral signatures, or to compare images taken at different times with different Sun and satellite positions and angles. By applying the corrections described here, we transform the satellite-derived reflectance into values that are (almost) directly comparable with those we would obtain with a field spectrometer. I wrote “almost” because we are not yet considering atmospheric correction, and also for two further assumptions we will make below. Also note that topographic correction just compensates for differences in measured reflectance due to terrain, as if they were lying flat and measured with a nadir-looking instruments, but cannot recreate values for areas behind obstructions.

We will need the calibrated images (radiance top-of-atmosphere, TOA, see previous posts); information about the sun elevation (Sun Zenith) and direction (Sun Azimuth), both from the MOD03 metadata file; information on the solar irradiance on each band, from each MOD02 hdf file; and two auxiliary grids, covering exactly the same region as the images, and containing the slope and aspect maps of the region of interest. The last two grids are derivable from a digital elevation model using standard GIS procedures (see, for example, GRASS).

The basic relation we will consider here is:

$\rho_i = \frac{\pi L_i^{TOA}}{E_i cos\beta cos\gamma}$

where $\rho_i$ is the corrected band reflectance, $L_i^{TOA}$ is the radiance at top-of-atmosphere, $E_i$ is the bandwise solar irradiance, $\gamma$ is the sensor viewing angle, or sensor zenith, and $cos \beta$ is derived as follows:

$cos \beta = cos \theta_n cos \theta_s + sin \theta_n sin \theta_s cos(\phi_n - \phi_s)$

where $\theta_n$ is the terrain slope, $\theta_s$ is the solar zenith, $\phi_n$ is the terrain aspect, and $\phi_s$ is the solar azimuth.

For the sake of this example, we assume that the materials on the surface are Lambertian, that is, they are perfectly diffuse isotropic reflectors, spreading the energy homogeneously in all directions. Admittedly, this is not true for many important natural surfaces, but works fairly OK as a first approximation, and allows us to compensate for illumination effects by simply dividing the observed radiance by the cosine of the sensor zenith angle. Think that more realistic assumptions about the surface materials invariably involve very complicated additional calculations which, more often than not, rapidly push the problem to the edge of tractability. Here we will take the pragmatic approach of proposing a “good-enough” solution.

In this post we will use the following MODIS band 2 image of the Southern Patagonia Icefield as an example:

The image covers approximately 395 km in the North-South direction, and 119 km in the East-West direction.

We will need to extract the information on solar zenith and azimuth, as well as sensor zenith from the MOD03 files. Note that, as MODIS images cover a large geographical area, these angles cannot be considered constant over a scene, as we will sometimes assume for Landsat, for example. Therefore it is important to extract these values pixel-wise, as grids from the metadata file. We begin by reading the image and figuring out what is offered there:

read_sds_attributes -sds='SolarZenith,SolarAzimuth,SensorZenith' MOD03....

and we get:

====================================================================== SDS : SolarZenith Attribute        Data Type           Value -----------------+-----------+--------------------------------------- units               STRING    degrees valid_range         INT16     0, 18000 _FillValue          INT16     -32767 scale_factor        FLOAT64   0.010000 ====================================================================== SDS : SolarAzimuth Attribute        Data Type           Value -----------------+-----------+--------------------------------------- units               STRING    degrees valid_range         INT16     -18000, 18000 _FillValue          INT16     -32767 scale_factor        FLOAT64   0.010000 ====================================================================== SDS : SensorZenith Attribute        Data Type           Value -----------------+-----------+--------------------------------------- units               STRING    degrees valid_range         INT16     0, 18000 _FillValue          INT16     -32767 scale_factor        FLOAT64   0.010000 

Note that values range from 0 to 18000 and, as the scale value of 0.01 indicates, are expressed in hundredths of degrees. This is important to consider when performing the calculations.

We then proceed to extract these grids by executing the usual swath2grid command, as explained in previous posts. Images with sensor and sun angles extracted from MOD3 are always 1 Km pixel size, and will need to be resampled for precessing with MOD02QM or MOD02HM. Here we have, from left to right: Solar Azimuth, Solar Zenith and Sensor Zenith for our region:

Note how values vary across the pictures from low (dark) to high (bright).

We then continue by deriving the terrain angles from a digital elevation model. We use an elevation model compiled using data from Shuttle Radar Topography Mission (SRTM). The typical no-data “holes” in SRTM data were filled by creating a new grid through interpolation of the original values (as described in my recent paper “Hypsometry and sensitivity of the mass balance to changes in equilibrium-line altitude: the case of the Southern Patagonia Icefield“). The elevation, and the derived aspect and slope grids look like this:

Before finishing, we must obtain the solar irradiance at each band ($E_i$ in Eq. 1). These values can be exprted using gdalinfo:

gdalinfo MOD02*

What we get is something like:

Solar Irradiance on RSB Detectors over pi=511.46, 511.46, 511.46, 511.46, 511.46, 511.46, 511.46, 511.588, 511.62, 511.588, 511.588, 511.588, 511.556, 511.556, 511.556, 511.492, 511.492, 511.556, 511.588, 511.588, 511.588, 511.62, 511.62, 511.62, 511.62, 511.62, 511.556, 511.492, 511.46, 511.365, 511.301, 511.206, 511.11, 510.983, 510.824, 510.378, 509.614, 509.614, 509.614, 509.614, 315.763, 315.763, 315.763, 315.763, 315.763, 315.763, 315.763, 315.795, 315.795, 315.795, 315.795, 315.795, 315.795, 315.763, 315.795, 315.795, 315.795, 315.795, 315.795, 315.795, 315.795, 315.795, 315.795, 315.795, 315.795, 315.795, 315.763, 315.763, 315.763, 315.732, 315.732, ...

which is a list of the solar irradiance value on each of the 40 sensors divided by pi (that is the inverse of the first part of the right hand side of Eq. 1, including the pi on the numerator). Now, in our simple example we will restrict ourselves to make a simple average of the values for each band, but we note that rigorous calibration requires detailed correction for each detector. I do this by parsing the “global attributes” text file with Perl and then averaging sequentially over 40 values using Perl Data Language.

We are now ready to combine these grids. As described in the previous post, I use the NetCDF calculator grdmath, that comes with GMT. The resultant topographically corrected surface reflectance image looks like this:

Note how the image looks “flat” now, in comparison with the uncorrected image. This is a Float32 bit image, with reflectance values from 0 to 1. Errors in the procedures and values might cause some pixels to have values larger than 1 or smaller than 0. These cases need to be analyzed with care. Steep slopes and ridges are particularly difficult subjects, specially when covered with snow, as in my example. Always have a second look at your data and procedures!

The procedures that we have reviewed here are appropriate to make a quick and dirty topographic and illumination correction of MODIS images. There are, nevertheless, a number of delicate topics that we have not covered and that must be seriously studied in detail when performing rigorous correction. First, as stated above, we have neglected to correct for atmospheric effects; second, we have blindly assumed a Lambertian surface, and third, we have neglected the individual response of each of the 40 detectors in MODIS, averaging their registered solar irradiance. We also did not quantify errors in reflectance. So, in conclusion, although our results look fairly good to the eye and the pixel values seem very reasonable, remember that these procedures are limited.

All in all, I hope readers of these posts will find them useful and will find themselves motivated to play around with MODIS and may be bring their analysis to more rigorous levels. Good luck and enjoy!

NOTE: originally published in my old blog “Tales of Ice and Stone”, on 2014-03-05