Newsletters - MATLAB Digest
Reading and Visualizing Data from the Earth Observing System (EOS)
(Part 3 of 4)
Reading MODIS Sea Surface Temperature Swath
With our color mapping planned, we turn to reading and visualizing our selected data sets, beginning with MODIS Sea-Surface Temperature (SST). Notice how we can use hdfinfo before calling hdfread. For this example, we also subset the SST swath, mainly to reduce the amount of data that we ultimately pass to surface.
mod28L2File = 'MOD28L2.A2001185.0830.003.2001308112641.hdf';
sstInfo = hdfinfo(mod28L2File,'eos');
sstSwath = sstInfo.Swath;
idx = {[1 1],[4 4],[]}; % Subset by 4 in each dimension
sst = hdfread(sstSwath,'Fields','sst','Index',idx);
After reading the scaled and shifted SST array into sst, we read the data flags and construct a land mask using hdfread and bitget.
commonFlags = hdfread(sstSwath,'Fields','common_flags','Index',idx);
land = (bitget(commonFlags,8) == 1);
Next, we read the shifting and scaling parameters (slope and intercept) and other metadata elements that are not strictly part of the HDF-EOS structure.
sstInfoHDF = hdfinfo(mod28L2File);
datasets = sstInfoHDF.Vgroup.Vgroup(2).SDS;
sstDataset = datasets(strmatch('sst',{datasets.Name},'exact'));
sstAttr = sstDataset.Attributes;
sstSlope = double(sstAttr(strmatch('Slope', {sstAttr.Name})).Value)
sstIntercept = double(sstAttr(strmatch('Intercept',{sstAttr.Name})).Value)
sstUnits = sstAttr(strmatch('Units',{sstAttr.Name})).Value
Finally, we read the MODIS geolocation data fields, converting the null value (-999) to NaN. Our use of variable idx with parameter 'Index' ensures that subsetting is consistent with the SST data itself.
mod03File = 'MOD03.A2001185.0830.003.2001305034045.hdf';
geoInfo = hdfinfo(mod03File,'eos');
sstLat = double(hdfread(geoInfo.Swath,'Fields','Latitude', 'Index',idx));
sstLon = double(hdfread(geoInfo.Swath,'Fields','Longitude','Index',idx));
sstLat(sstLat == -999) = NaN;
sstLon(sstLon == -999) = NaN;
Mapping SST Values to the Colormap
The SST values are encoded in the HDF-EOS data file as unsigned 16-bit integers. (Array sst is MATLAB class uint16.) To convert to physical temperature units (degrees Celsius), we need to multiply by sstSlope, then add sstIntercept. However, instead of converting to an intermediate, double-valued array in degrees Celsius, we scale and shift the encoded SST values to correspond directly to the SST section of the colormap, in preparation for visualizing the SST product via direct color mapping. To do this, we determine coefficients a(1) and a(2) such that a(1) * sst + a(2) maps the encoded SST values directly onto the SST section of the colormap.
% Find the encoded SST limits that correspond to the temperature limits
% in sstDataLim.
sstCodelim = (sstDataLim - sstIntercept)/sstSlope;
% Solve a 2-by-2 linear system:
% a(1) * sstCodelim + a(2) = sstCmapLim - 1
a = (sstCmapLim - 1)/[sstCodelim; [1 1]];
We subtract 1 from sstCmapLim because we are going to use class uint8 to store our remapped SST values. For convenience, MATLAB automatically translates the range [0 255] to [1 256] when using indexed color with 8-bit imagery.
Finally, we apply the coefficients in array a, clip to stay within the appropriate part of the colormap, handle special values, and convert to uint8.
sstMapped = a(1) * double(sst) + a(2); % Scale and shift
sstMapped(sstMapped < sstCmapLim(1)-1) = sstCmapLim(1)-1; % Clip bottom
sstMapped(sstMapped > sstCmapLim(2)-1) = sstCmapLim(2)-1; % Clip top
sstMapped(sst == 0) = 0; % Remap special value (displays as white)
sstMapped(land) = 1; % Apply land mask (displays as dark gray)
sstMapped = uint8(round(sstMapped)); % Save storage with uint8
Visualizing the MODIS SST Swath
Having remapped the SST data, we display it as an indexed-color image:
figure('Position',[100 100 size(sstMapped,2) size(sstMapped,1)])
colormap(cmap)
image(sstMapped)
set(gca,'Position',[0 0 1 1],'DataAspectRatio',[1 1 1]);
 |
| Click on image to see enlarged view. |
The horizontal bar in the lower right is an artifact due to a sequence of records with partial dropouts or other quality problems.
Geometrically, the swath is simply a time-ordered sequence of records from the MODIS sensor, and hence includes several types of spatial distortion. We illustrate this by contouring the latitude and longitude geolocation arrays and displaying the contours over the image. We use a contour interval of 5 degrees:
interval = 5;
latContourValues = interval ...
* (ceil(min(sstLat(:))/interval) : floor(max(sstLat(:))/interval));
lonContourValues = interval ...
* (ceil(min(sstLon(:))/interval) : floor(max(sstLon(:))/interval));
hold on
[c,h] = contour(sstLat,latContourValues,'k'); clabel(c,h)
[c,h] = contour(sstLon,lonContourValues,'k'); clabel(c,h)
hold off
 |
| Click on image to see enlarged view. |
We use the geolocation fields as x and y inputs to surface to redisplay the swath correctly positioned with respect to latitude-longitude coordinates. We also use colorscale to add a color scale bar indicating degrees Celsius:
figure
colormap(cmap)
axes('DataAspectRatio',[1 1 1],...
'XTick',-180:5:180,'YTick',-90:5:90,...
'XLim',[20 50],'YLim',[12 37],...
'XGrid','on','YGrid','on')
xlabel('Longitude, degrees')
ylabel('Latitude, degrees')
title('Geolocated SST swath, July 4, 2001')
surface(sstLon, sstLat, double(sstMapped) + 1,...
'Linestyle','none','CDataMapping','Direct')
% Add color scale
pos = get(gca,'Position');
set(gca,'Position',pos + [-0.06 0 0 0])
colorscale(sstCmapLim, sstDataLim, 5, 'vert',...
'Position',[0.92 pos(2) 0.025 pos(4)])
xlabel('SST')
ylabel('degrees Celsius')
 |
| Click on image to see enlarged view. |
