Sensors operating in the visible and mid-infrared portions of the electromagnetic spectrum measure the fraction of the incoming solar radiation reflected at the surface. If the sensors are mounted on board satellites in orbit, the radiation coming up to the sensor does not entirely arises from the reflection at a particular location on Earth, but also partly due to the interaction of light with the atmosphere. Apart from the obvious case of a cloudy sky, absorption and dispersion of light by the atmosphere can be very important, even in conditions that we could consider clear sky.

The basic types of scattering are Rayleigh and Mie. Rayleigh scattering arises from the interaction between the light and the air molecules. it affects more the shorter wavelengths and is responsible for the blue color of the sky. Mie scattering arises from the interaction of light with particles, more importantly water drops, dust and pollution. Their combined effect on a satellite image is that, apart from the radiation reflected at the surface, we also get radiation coming from Sun light scattered in the atmosphere along the view path (path irradiance), and from Sun light scattered in the atmosphere in the direction of the surface element and reflected to the sensor (sky irradiance). Eventually we also get light scattered from the neighboring surface (neighbor irradiance).

Depending on our purpose these effects might be negligible and we can simply disregard them, as when measuring the extent of a feature. In other cases they may severely degrade the quality of an image, for example when we need to use the value of the surface reflectance. We might then need to apply atmospheric correction. In essence, this is about modeling the amount of radiation arising from the above described processes and subtracting them from the measured values. Atmospheric correction models, such as 6s, used in GRASS GIS, are based on radiative transfer theory and depend on information about current atmospheric conditions at the moment of image acquisition, be them measured (rarely) or assumed (most often). We here show a simple application of the haze correction method, a simple yet effective method that relies on image values and simple assumptions about the effects of the atmosphere for visible and near infrared Landsat bands. In the tradition of this blog, we focus on the fundamental aspects, leaving sophistication aside. As an example, we use a subset of a Landsat 8 image of the Bay of Naples, Italy, subject of previous posts.

The central assumption of the method is that all values in an image are shifted upwards due to the additive effects of scattering, and that the effects are uniformly distributed across the region covered by the image. Thus, if we could find the value of this upward shift in each band, we could simply subtract it from each pixel and correct the image. To find these values we plot histograms of the bands and inspect them to find the minimum significant values in each band. We begin with the reflectance in bands 2 (blue), 3 (green), 4 (red) and 5 (near infrared):

Histograms of reflectance in bands 2, 3, 4 and 5, of the Landsat 8 image path/row 189/32.
Histograms of reflectance in bands 2, 3, 4 and 5, of the Landsat 8 image path/row 189/32.

We clearly see that band 2, blue, has the largest shift in the histogram, the signature of Rayleigh scattering. Detailed inspection of the data for each histogram enables us to figure out the thresholds for the significant values:

band 2: 0.15
band 3: 0.11
band 4: 0.07
band 5: 0.04

We then subtract these threshold values from each band. We do this using grdmath, the NetCDF calculator of the Generic Mapping Tools (GMT) toolbox:

grdmath vesuvius_b2_refTOA.grd 0.15 SUB = vesuvius_b2_refTOA_atcorr.grd
grdmath vesuvius_b3_refTOA.grd 0.11 SUB = vesuvius_b3_refTOA_atcorr.grd
grdmath vesuvius_b4_refTOA.grd 0.07 SUB = vesuvius_b4_refTOA_atcorr.grd
grdmath vesuvius_b5_refTOA.grd 0.04 SUB = vesuvius_b5_refTOA_atcorr.grd

Note that, to do this, we first have to convert our GeoTIFF files to NetCDF using the GDAL utility gdal_translate. See previous posts. A comparison of our results, in band 2, with the uncorrected image readily shows the effects of “haze removal”:

Band 2 (blue), before (left) and after (right) removing the 'haze'
Band 2 (blue), before (left) and after (right) removing the ‘haze’

A comparison of false color composites, using bands 5, 4 and 3 (RGB), shows the same but in color:

False color composite (bands 5, 4 and 3 RGB), before (left) and after (right) haze removal.
False color composite (bands 5, 4 and 3 RGB), before (left) and after (right) haze removal.

Finally, we might calculate the Normalized Difference Vegetation Index (NDVI):

NDVI = \frac{\rho_{ir} - \rho_{r}}{\rho_{ir} + \rho_{r}}

using our corrected images, and compare with the results with the uncorrected ones:

Normalized Difference Vegetation Index (NDVI) calculated before (left) and after (right) haze removal.
Normalized Difference Vegetation Index (NDVI) calculated before (left) and after (right) haze removal.

Although there is no big visual difference, the image to the right is arguably better because the range of values in the corrected bands is larger, thus allowing us to reveal vegetation patterns more clearly.

This method is closely related to the Dark Object Subtraction (DOS) method, which is based on the identification of dark surfaces, clean lakes or asphalt, and assuming that their reflectances are only due to atmospheric scattering, subtracting their values from the image. The difference here is that we do not have to rely on the existence of such surfaces within our region. In a sense, haze removal by histogram minimum is technically equivalent to a linear stretch, where we cut a fraction of the lowest values, except that in this case we have a reasonable geophysical principle to rely on. Although not perfect, the haze removal by histogram minimum is a robust and simple method that is very reliable as a first approximation. In many cases it might be the only atmospheric correction we need.