Newsletters - MATLAB Digest

Reading and Visualizing Data from the Earth Observing System (EOS)

(Part 2 of 4)

Data Selection

We have chosen to work with two remote sensing data sets that provide information about the global environment in early July 2001. Each data set is an instance of a standard NASA data product.

The first is an HDF-EOS swath representing sea surface temperature (SST) derived from the Moderate Resolution Imaging Spectrometer (MODIS), one of the instruments on the Earth Observing System's Terra satellite. (A second MODIS instrument went into orbit this spring on board the next satellite in the EOS series, Aqua.) The MODIS swath extends over the area covered during a five-minute period on July 4, 2001. It covers parts of the eastern Mediterranean Sea and the northern Red Sea, at one-kilometer resolution.

The second data set is a one-degree global grid representing the state of the world's vegetation via a quantity known as a Normalized Difference Vegetation Index (NDVI). Each one-degree bin averages results from the Advanced Very High Resolution Radiometer (AVHRR) on the National Oceanic and Atmospheric Administration (NOAA) Polar Orbiter satellite NOAA-14 for the period of July 1-10,2001. This is an instance of a much older data product, and is distributed in HDF rather than HDF-EOS.

Both the SST and NDVI data products take advantage of the multispectral nature of the sensors from which they were derived. The SST product gauges the temperature of waters near the ocean surface using thermal infrared bands bracketing a wavelength of 11 microns (MODIS bands 31 and 32). The principle underlying NDVI is the strong reflectance of healthy, chlorophyll-based vegetation at near infrared wavelengths and its relatively weak reflectance in the visible red. NDVI is simply defined as the difference between these reflectances normalized by their sum. This yields a dimensionless quantity ranging in theory from -1 to 1 (but in practice seldom falling below -0.25.

NASA's Distributed Active Archive Centers (DAACs) are the primary sources for Earth Observing System data products. The Goddard Space Flight Center (GSFC) Earth Sciences DAAC emphasizes products related to the earth's atmosphere and global biosphere. We ordered all the data sets used in this article from this DAAC. See daac.gsfc.nasa.gov for rich sources of information on SST, NDVI, and many other data products.

We obtained the following data files:

Filename	File Type	Contents
MOD28L2.A2001185.0830.003.2001308112641.hdf	HDF-EOS	MODIS SST swath and metadata
MOD03.A2001185.0830.003.2001305034045.hdf	HDF-EOS	MODIS geolocation fields
PAL_CLIMATE_JUL_01-10_2001.HDF	HDF	AVHRR NDVI grid and metadata

The first file contains the MODIS SST swath (and metadata such as quality masks, text descriptions, etc.). It is complemented by the second file, which contains two geolocation fields: an array of latitudes and an array of longitudes that combine to give the geographic location of every element in the swath. The third file contains the AVHRR NDVI product and its metadata. This data product was processed so that elements are regularly positioned in latitude and longitude; there is no need for additional geolocation data.

Exploring with HDFTOOL

We can explore our data files with the HDF Import Tool:

hdftool('MOD28L2.A2001185.0830.003.2001308112641.hdf')
hdftool('MOD03.A2001185.0830.003.2001305034045.hdf')
hdftool('PAL_CLIMATE_JUL_01-10_2001.HDF')

The first call starts the tool and the next two calls detect an open hdftool and add more data files to it.

Click on image to see enlarged view.

Looking at the HDF Import Tool window, the drop-down View menu says "File as HDF and HDF-EOS", and the two HDF-EOS files (MOD28L2... and MOD03...) appear in both of the panels on the left side of the window. The HDF file (PAL_CLIMATE...) appears only in the "File viewed as HDF" window.

In the "File viewed as HDF-EOS" panel, both HDF-EOS data files have a Geolocation Fields object, but the first one is empty. This confirms our need for the MOD03 product to supplement the MOD28L2 product.

For the MOD28L2 product, the swath data fields appear in both the HDF and HDF-EOS views. The data fields of concern are "sst" (with the 11-micron SST) and "common_flags" (with a bit to mask out land areas). The "sst" field is highlighted under the HDF view, and the corresponding metadata are displayed in hdftool in the upper right panel. In addition to the swath dimensions, we see that the sea surface temperature values are in units of degrees Celsius, but are encoded as unsigned 16-bit integers. A simple linear formula for decoding (converting the value in the file to a value in degrees Celsius) is given, along with the required slope and intercept values.

Note: In the MOD28L2 file, the slope and intercept metadata are stored in additional HDF data sets within the HDF-EOS file. This is permissible, but such fields are not visible from the HDF-EOS view (nor when using only the official HDF-EOS interface libraries). Therefore, it is important that we use the HDF view in addition to the HDF-EOS view.

When bringing the data sets into MATLAB, we use calls to hdfinfo and hdfread. This way we illustrate these additional functions and we keep our narrative more streamlined. However, it is easy to import data interactively; that is the purpose of the large panel in the lower-right corner of the HDF Import Tool window.

Color Scales and Sharing a Colormap

One of the key goals here is to display the SST and NDVI data products in the same figure, in their correct geographic relationship. Indexed color is a natural choice for visualizing both data sets, but we have a problem: The color scales traditionally used for SST and NDVI are very different and MATLAB permits only one colormap per figure. One solution (the one we'll adopt) is to create a colormap with separate portions for SST and NDVI, and then manage which part of the color map is used by scaling and shifting the data sets ourselves and setting 'CDataMapping' to 'direct' in calls to image and surface. (A possible alternative would be to apply the NDVI colormap to create an RGB image, and use indexed color only for the SST data.)

With help from the new Colormap Editor in MATLAB, we created a 256-element colormap (a 256-by-3 array of class double) to support both SST and NDVI and saved it in a MAT file ('colormaps.mat'). We will use this colormap for every figure in this article. The first half of the colormap is dedicated to SST and the second half is for NDVI data. In general, a shared colormap such as this can be much longer, can be unevenly divided among data sets, and may support more than just two data sets. One advantage of the choice we made is that we can economize on storage by fitting the scaled and shifted data into unsigned 8-bit integers (MATLAB class uint8).

Of course, remapping the data to uint8 can cause both round off and loss of dynamic range, especially in the case of the SST data that originates as unsigned 16-bit integers. This effect is acceptable for visualization, which is our purpose. However, we would have to return to the original data (and possibly convert to physical units and MATLAB class double) before performing any precise quantitative analyses.

We can easily load our colormap into MATLAB and display it using the colorbar function:

load colormaps
cmap = [sstCmap; ndviCmap]; % We saved our colormap in two equal parts
figure
colormap(cmap)
set(axes,'Visible','off')
h = colorbar('horiz');
set(get(h,'Title'),'String','Standard Colorbar')

Click on image to see enlarged view.

It is also easy to examine the colormap with the Colormap Editor:

colormapeditor

Click on image to see enlarged view

Looking at the Colormap Editor window shown here, notice that the cursor is at 130, a pure gray (see the RGB and HSV readouts in the middle of the screen) that we use for null data areas in the NDVI product. To its left is a deep blue used for water areas. A continuous color scale for NDVI extends to the right from 131 onward. The upper end of the SST scale appears on the left, ending at 128. Special values at 1 for null or missing data (white) and 2 for land (dark gray) appear outside the window on the far left end of the scale.

Our next step is to define four 1-by-2 arrays to be used throughout. They assign ranges within the colormap to SST and NDVI, along with the corresponding ranges in physical units:

sstCmapLim = [3 128];
sstDataLim = [0 35];
ndviCmapLim = [ 131 256];
ndviDataLim = [-0.25 1.00];

The array sstCmapLim designates the continous interval on the colormap to be used to visualize actual temperature values. sstDataLim designates the corresponding values at these limits (in this case, degrees Celsius). ndviCmapLim and ndviDataLim have similar meanings, except that the physical units are dimensionless. Each pair of arrays provides the first two parameters for our colorscale function, which takes the following syntaxes:

H = COLORSCALE(CMAPLIM,DATALIM,DATAINC,ORIENTATION,'Position',RECT)
H = COLORSCALE(...,PROP1,VAL1,PROP2,VAL2,...)

The data increment (DATAINC) controls the placement of ticks on the color scale axes. The orientation should be either 'horiz' or 'vert'. Any property-value pair that can be used with axes is accepted, and a 'Position' property is required. Unlike colorbar, colorscale requires that you designate a position within the current (or parent) figure.

The next step is to add color scales for SST and NDVI to the current figure using colorscale and our CmapLim and DataLim arrays. Notice how we use 'Position' to place the color scale in the figure with respect to an existing axes (in this case, the one created with colorbar). We use this programming pattern in all our examples.

pos = get(h,'Position'); % Position color scales relative to the color bar
colorscale(sstCmapLim, sstDataLim,5,'horiz','Position',pos + [0 0.70 0 0])
title('SST Colorscale')
xlabel('degrees Celsius')
colorscale(ndviCmapLim,ndviDataLim,0.2,'horiz','Position',pos + [0 0.35 0 0])
title('NDVI Colorscale')

Click on image to see enlarged view

Part 3