In the third part of this series we will calibrate the images to convert the pixel values to spectral radiance at top of atmosphere.

The spectral radiance is a physical quantity expressing the amount of energy being radiated (in our case, back to space), and is measured in Watts (Joules per second, energy flux) per square meter (area) per micrometer (wavelength) per steradian (solid angle).

The values that we will obtain here will not yet be corrected for illumination, topographic or atmospheric effects, and therefore correspond to what is measured by the satellite at the altitude of its orbit, as opposed to values on the ground that you can measure with a field spectrometer. The values we will obtain here are thus called “radiance at top-of-atmosphere”, usually abbreviated “TOA”.

The conversion is based on a linear relationship between data counts (pixel values) and the radiance measured by the sensor (for the gory details have a look at the manuals in NASA’s MODIS Characterization Support Team). The basic relationship is:

radiance = (pixel_value - offset) * scale

and has to be performed band-wise. We begin by extracting the scales and offset from bands 3-7, that is the 500 m bands (remind that we use the bands resampled to 1 km). We use the program read_sds_attributes that comes with the MODAT suite.

read_sds_attributes -sds='EV_500_Aggr1km_RefSB' MOD021KM.A2013187.1035.005.2013187194950.hdf

and we get:

======================================================================
SDS : EV_500_Aggr1km_RefSB
Attribute Data Type Value
-----------------+-----------+---------------------------------------
long_name STRING Earth View 500M Aggregated 1km
Reflective Solar Bands Scaled
Integers
units STRING none
valid_range UINT16 0, 32767
_FillValue UINT16 65535
band_names STRING 3,4,5,6,7
radiance_scales FLOAT32 0.035233, 0.024804, 0.005951, 0.002677, 0.000823
radiance_offsets FLOAT32 -0.000000, -0.000000, -0.000000, -0.000000, -0.000000
radiance_units STRING Watts/m^2/micrometer/steradian

reflectance_scales FLOAT32 0.000055, 0.000043, 0.000041, 0.000036, 0.000030
reflectance_offsets FLOAT32 -0.000000, -0.000000, -0.000000, -0.000000, -0.000000
reflectance_units STRING none
corrected_counts_scalesFLOAT32 0.124973, 0.124973, 0.124973, 0.124973, 0.124973
corrected_counts_offsetsFLOAT32 -0.000000, -0.000000, -0.000000, -0.000000, -0.000000
corrected_counts_unitsSTRING counts
Processing done !

The output expresses the scales (reflectance_scales) and offsets (reflectance_offsets) for each of the five resampled bands. The same is then done for bands 1 and 2, the 250 m bands:

read_sds_attributes -sds='EV_250_Aggr1km_RefSB' MOD021KM.A2013187.1035.005.2013187194950.hdf

======================================================================
SDS : EV_250_Aggr1km_RefSB
Attribute Data Type Value
-----------------+-----------+---------------------------------------
long_name STRING Earth View 250M Aggregated 1km
Reflective Solar Bands Scaled
Integers
units STRING none
valid_range UINT16 0, 32767
_FillValue UINT16 65535
band_names STRING 1,2
radiance_scales FLOAT32 0.026587, 0.009931
radiance_offsets FLOAT32 -0.000000, -0.000000
radiance_units STRING Watts/m^2/micrometer/steradian

reflectance_scales FLOAT32 0.000054, 0.000033
reflectance_offsets FLOAT32 -0.000000, -0.000000
reflectance_units STRING none
corrected_counts_scalesFLOAT32 0.124973, 0.124973
corrected_counts_offsetsFLOAT32 -0.000000, -0.000000
corrected_counts_unitsSTRING counts
Processing done !

We can then proceed to extract these values and recalculate the images. I find convenient to do these procedures via Perl scripts, by executing read_sds_attributes command, scanning the output text file, and then extracting the values using regular expressions. I then perform the calculations by converting the GeoTIFFs to NetCDF grids and calling the grdmath tool that comes with Generic Mapping Tools:

grdmath file_raw.grd offset SUB scale MUL = file_corrected.grd

* Note: grdmath is a Reverse Polish Notation grid calculator.

and lastly call the gdal_translate utility program from GDAL to produce a new GeoTIFF:

gdal_translate file.grd file.tif

Since the calibration is a linear conversion, the results look pretty much the same as before on the screen. The only difference is that our pixel values are now geophysical quantities, as seen in the following example:

hps

The images depict the uncorrected pixel counts (left) and radiance at top-of-atmosphere (right) for the same region, the southernmost portion of the Southern Patagonia Icefield, as seen in MODIS band 2 (841 – 876 nm). The mouse pointer is located over (roughly) the same place over the ablation area of Tyndall Glacier. The pixel values can be read on the lowermost part of the image.

We have here reviewed how to derive TOA radiance values from uncorrected pixel counts in MODIS L1B images. We did not discuss how to evaluate errors, or take into account uncertainties inherent to the sensors. May be we can handle that in a future post. In the next and last post, we will see how to perform illumination and topographic correction to achieve a product that is suitable for multitemporal comparison and analysis.

NOTE: originally published in my old blog Tales of Ice and Stone, on 2014-02-20

Advertisements