Newsletters - MATLAB Digest
Reading and Visualizing Data from the Earth Observing System (EOS)
(Part 4 of 4)
Reading AVHRR Global NDVI Grid
From exploratory work with hdftool and hdfinfo, we know that the NDVI data is stored in 'Data-Set-2' within the NDVI file, and we read it directly into a MATLAB array using hdfread:
ndviFile = 'PAL_CLIMATE_JUL_01-10_2001.HDF';
ndvi = hdfread(ndviFile,'Data-Set-2');
We use hdfinfo to get the offset and scale factor attributes that relate the integers in the file to absolute NDVI values:
ndviInfo = hdfinfo(ndviFile);
ndviAttr = ndviInfo.SDS.Attributes;
ndviAddOffset = double(ndviAttr(...
strmatch('add_offset', {ndviAttr.Name},'exact')).Value)
ndviScaleFactor = double(ndviAttr(...
strmatch('scale_factor',{ndviAttr.Name},'exact')).Value)
Mapping NDVI Values to the Colormap
As with the SST data, the NDVI values are encoded in the HDF file as integers (except that they are unsigned 8-bit instead of unsigned 16-bit). Therefore we take the same approach to remapping the data, taking into account that the linear relationship between the encoded NDVI values and actual NDVI levels is expressed in terms of an offset and scale factor rather than an intercept and slope. To convert an encoded data value to absolute NDVI units, we need to subtract the offset, then multiply by the scale factor.
% Find the encoded NDVI limits that correspond to the absolute NDVI limits
% in ndviDataLim.
ndviCodelim = ndviDataLim/ndviScaleFactor + ndviAddOffset;
% Solve a 2-by-2 linear system:
% b(1) * ndviCodelim + b(2) = ndviCmapLim - 1
b = (ndviCmapLim - 1)/[ndviCodelim; [1 1]];
Here again, we subtract 1 in order to store the remapped values as uint8.
Finally, as before, we apply the coefficients in b, clip to stay within the appropriate part of the colormap, handle special values, and convert to uint8.
ndviMapped = b(1) * double(ndvi) + b(2); % Scale and shift
ndviMapped(ndviMapped < ndviCmapLim(1)-1) = ndviCmapLim(1)-1; % Clip bottom
ndviMapped(ndviMapped > ndviCmapLim(2)-1) = ndviCmapLim(2)-1; % Clip top
ndviMapped(ndvi==0) = 130-1; % Missing data (displays as gray)
ndviMapped(ndvi==1) = 129-1; % Water (displays as blue)
ndviMapped = uint8(round(ndviMapped)); % Save storage with uint8
Visualizing the NDVI Grid
Having remapped the NDVI grid to match the colormap, we display it as an image. We use global latitude-longitude axes and include a color scale:
figure
colormap(cmap)
image(ndviMapped,'XData',[-179.5 179.5],'YData',[89.5 -89.5])
set(gca,'YDir','normal','DataAspectRatio',[1 1 1],...
'XTick',[-180 : 30 : 180],'YTick',[-90 : 30 : 90],...
'XLim',[-180 180],'YLim',[-90 90])
xlabel('Longitude, degrees')
ylabel('Latitude, degrees')
title('NDVI from AVHRR, July 1-10, 2001')
% Add color scale
pos = get(gca,'Position');
set(gca,'Position',pos + [0 0.1 0 0])
colorscale( ndviCmapLim, ndviDataLim, 0.2, 'horiz','Position',[pos(1:3) 0.05])
xlabel('NDVI')
 |
| Click on image to see enlarged view. |
Combined SST and NDVI Display
Next we display the NDVI and SST data together. Our latitude-longitude window includes the entire swath and the general area surrounding it. We add a color scale for each data set.
figure('Renderer','zbuffer')
colormap(cmap)
image(ndviMapped,'XData',[-179.5 179.5],'YData',[89.5 -89.5])
set(gca,'YDir','normal','DataAspectRatio',[1 1 1],...
'XTick',-180:5:180,'YTick',-90:5:90,...
% Ticks every 5 degrees
'XLim',[15 60],'YLim',[0 45]) % Axes limits in degrees
xlabel('Longitude, degrees')
ylabel('Latitude, degrees')
title('SST swath plotted over one-degree NDVI grid')
surface(sstLon, sstLat, ones(size(sstLon)), double(sstMapped) + 1,...
'Linestyle','none','CDataMapping','Direct');
% Add color scales
pos = get(gca,'Position');
set(gca,'Position',pos + [-0.12 0 0 0])
colorscale(ndviCmapLim, ndviDataLim, 0.2, 'vert',...
'Position',[0.9 pos(2) 0.025 pos(4)])
xlabel('NDVI')
colorscale(sstCmapLim, sstDataLim, 5, 'vert',...
'Position',[0.8 pos(2) 0.025 pos(4)])
xlabel('SST')
ylabel('degrees Celsius')
 |
| Click on image to see enlarged view. |
Map Projection of Combined SST and NDVI
With the Mapping Toolbox, we can move beyond displaying data over straight latitude-longitude coordinates---and the distortions that causes---to a wide range of projected map coordinate systems. For example, we can choose a projection to eliminate local shape distortions or to ensure that features of equal area on the globe have equal area on the map. Here we use the Mapping Toolbox commands axesm, meshm, and surfm to display the NDVI and SST data in a Mollweide (pseudocylindrical) equal-area projection.
First we set up our figure, title, colormap, and map axes, choosing a Mollweide projection with Greenwich as its central meridian.
figure('Color','k');
title({'SST from MODIS and NDVI from AVHRR, July 2001'},'Color','w')
colormap(cmap);
ax = axesm('MapProjection','mollweid','Origin',[0 0 0]);
set(ax,'Color',[0 0 0])
We display the global NDVI grid using meshm. We give this Mapping Toolbox function a "regular matrix map," with columns running from south to north and positioned regularly in latitude and longitude.
ndviNorthwestCorner = [89.5 -179.5]; % Lat/lon of upper-left-most cell center ndviCellSize = 1; % Size of grid cells in degrees
hmesh = meshm(flipud(double(ndviMapped) + 1),... % Converting uint8 image
[ndviCellSize ndviNorthwestCorner],... % The 'legend'
'CDataMapping','direct'); % Go directly into the colormap
We display the MODIS SST swath over the global NDVI grid using surfm, which takes the geolocation arrays as its first two arguments. The Mapping Toolbox considers the swath to be a "generalized matrix map."
hsurf = surfm(sstLat,sstLon,double(sstMapped) + 1,...
'CDataMapping','direct');
Finally, we add a color scale for each data set.
set(ax,'Position',get(ax,'Position') + [0 0.12 0 0])
colorscale(sstCmapLim,sstDataLim,5,'horiz',...
'Position',[0.1 0.10, 0.8, 0.03],'XColor','w','YColor','w')
title('Sea Surface Temperature','Color','w')
xlabel('degrees Celsius')
colorscale(ndviCmapLim,ndviDataLim,0.2,'horiz',...
'Position',[0.1 0.28, 0.8, 0.03],'XColor','w','YColor','w')
title('Normalized Difference Vegetation Index','Color','w')
 |
| Click on image to see enlarged view. |
Conclusions
Support for the HDF and HDF-EOS formats, plus graphics and visualization capabilities, make MATLAB an attractive environment to explore and display earth remote sensing data from Terra/MODIS, AVHRR, and other Earth Observing System satellites and sensors. The new HDF Import Tool and Colormap Editor GUIs make MATLAB 6.5 particularly useful. The function colorscale provides additional options for combining multiple data sets where each has its own color scale. The Mapping Toolbox enables displays in a wide range of projected map coordinate systems. Taken together, these tools provide a flexible, highly effective system for importing and visualizing EOS data products.
Files required to run this demonstration can be downloaded at:
www.mathworks.com/matlabcentral/fileexchange/Files.jsp?fileId=2611.
To learn more about products mentioned in this article, visit:
MATLAB
Mapping Toolbox
