How to Read and Visualize HDF-EOS5 data Using MATLAB

MATLAB provides APIs for accessing HDF5 data including HDF-EOS5. To use these APIs, user must be familiar with the HDF5 library APIs. This page presents a few examples on how to read and visualize HDF-EOS5 data via MATLAB.

Grid

HDF-EOS5 grid data is special in that it does not have to store longitude and latitude values. Instead of storing those values, HDF-EOS5 grid may keep only several parameters that can generate the entire longitude and latitude values. This can save some disk space, but this also makes it difficult to retrieve those values unless HDF-EOS5 API is used.

The degree of difficulty of obtaining longitude and latitude values depends on the projection method used in HDF-EOS5 grid. The simplest projection is the geographic projection. We will cover this case using a real OMI level 3 Ozone external data file from GES DISC. You can download the file here.

To access a data field in a grid data, the HDF-EOS5 file should be opened using H5F.open, first. Using the descriptor returned by that API, the grid should be opened. In this example, we will access a grid dataset called ColumnAmountO3 in an HDF-EOS5 file named OMI-Aura_L3-OMTO3e_2017m0105_v003-2017m0203t091906.he5.

Figure 1 Accessing a grid

      FILE_NAME='OMI-Aura_L3-OMTO3e_2017m0105_v003-2017m0203t091906.he5';

      DATAFIELD_NAME = '/HDFEOS/GRIDS/OMI Column Amount O3/Data Fields/ColumnAmountO3';

      file_id = H5F.open(FILE_NAME, 'H5F_ACC_RDONLY', 'H5P_DEFAULT');

      data_id = H5D.open(file_id, DATAFIELD_NAME);

H5D.open is the function which MATLAB provides to read data from an HDF-EOS5 grid file.

Then, dimensions, datasets, and attributes in this grid can be accessed.

Figure 2 Reading a data field in a grid

      data_space = H5D.get_space(data_id);

      [data_numdims data_dims data_maxdims] = H5S.get_simple_extent_dims(data_space);

      data = H5D.read(data_id,'H5T_NATIVE_DOUBLE', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT');

      ATTRIBUTE = '_FillValue';

      attr_id = H5A.open_name (data_id, ATTRIBUTE);

      fillvalue=H5A.read (attr_id, 'H5T_NATIVE_DOUBLE');

The data in a grid can be read out by using H5D.read. Note that data_id, is the descriptor returned by H5D.open. data is the buffer where the data is stored. Unlike using C APIs, users do not need to allocate memory for this buffer; memory will be automatically allocated. Also, users do not need to specify the type of data because MATLAB is a dynamically typed language.

After finishing reading all fields, one needs to close the grid using H5D.close. Similarly, if accessing the grid file is done, one needs to close the file using H5F.close.

Figure 3 Closing the grid and the file

      H5S.close (data_space);
      
      H5A.close(attr_id);

      H5D.close(data_id);

      H5F.close(file_id);

data_id is what H5D.open returned, and file_id is what H5F.open returned.

A few more steps are required to visualize the data field we just read. To generate a more meaningful plot, the data field needs to be associated with coordinates. Users need to provide longitude and latitude values. Note that the following explanation only applies to the geographic projection.

The geographic projection, which is used by OMI Column Amount O3 grid, maps meridians to equally spaced vertical straight lines, and circles of latitude to evenly spread horizontal straight lines ^[1]. This implies that we can precisely interpolate all longitude and latitude values if we know the followings:

leftX, rightX: the range of longitude
numX: the number of points in X dimension
offsetX: offset related to Pixel Registration
upperY, lowerY: the range of latitude
numY: the number of points in Y dimension
offsetY: offset related to Pixel Registration

From the range and the number of points, we can get the space between two adjacent points, which can be called the scale. Then, we can say i-th longitude value is (i + offsetX) * scaleX + leftX. offsetX is given by Pixel Registration. This is out of scope of this page, the value of offsetX and offsetY can be 0, 0.5 or 1. For more information, please refer to the HDF-EOS5 Reference ^[2].

We need to retrieve those four types of information. We will use HDFView to get the values for all the required variables. Using HDFView, you can open StructMetadata.0 string dataset under HDFEOS INFORMATION group. The string data contains XDim, YDim, UpperLeftPointMtrs, LowerRightMtrs, PixelRegistration, and GridOrigin information.

Figure 4 Values of the variables needed to plot the grid data

The usage of the above values for calculating latitude and longitude values can be seen in Figure 5.

Figure 5 Setting lat/lon information about a grid structure

      offsetY = 0.5;

      offsetX = 0.5;

      scaleX = 360/data_dims(2);

      scaleY = 180/data_dims(1);

      for i = 0:(data_dims(2)-1)

          lon_value(i+1) = (i+offsetX)*(scaleX) + (-180);

      end

      for j = 0:(data_dims(1)-1)

          lat_value(j+1) = (j+offsetY)*(scaleY) - 90;

      end

data_dims(2) (1440) and data_dims(1) (720) are the values for numX and numY, respectively. Also, -180 and -90 represents coordinates of upper left corner (UpperLeftPointMtrs) and lower right corner (LowerRightMtrs) respectively in DMS format^[3].

Finally, one can use the contour function to draw a plot using data, lon and lat calculated above. Since Many options are provided for using the surfm function, users may need to refer to MATLAB's detailed document for more information about this function.

Figure 6 Visualizing a data field

      surfm(lon, lat, data);

      coast = load('coast.mat');

      plotm(coast.lat,coast.long,'k');

You can see the complete code from here. Figure 7 shows the result.

Contour plot for a data field, RrLandRain

Figure 7 Plot for a data field, O3

Reference

Equirectangular projection. Wikipedia. http://en.wikipedia.org/wiki/Geographic_projection
HE5_GDdefpixreg page 2-199 HDF-EOS Interface Based on HDF5, Volume 2: Function Reference Guide Raytheon Company
6.5.4 Projection Parameters HDF-EOS Interface Based on HDF5, Volume 1: Overview and Examples Raytheon Company

Last modified: 11/11/2020