Working with MODIS L1B from scratch. 4. topographic and illumination correction

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: 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: