|
|
|
MATLAB
provides APIs for accessing several
HDF4
application programming interfaces (APIs) including
HDF-EOS2
. In order to use these APIs user must be familiar with the HDF4 library. This page presents a few examples on how to read and visualize HDF-EOS2 data via MATLAB.
MATLAB provides functions to read HDF-EOS grid (GD), point (PT), and swath (SW) data. Since each of them is equivalent to an HDF-EOS2 C API, those who are familiar with the HDF-EOS2 library can easily learn how to handle HDF-EOS2 files in MATLAB.
HDF-EOS grid data is special in that it does not store longitude and latitude values. Instead of storing those values, HDF-EOS grid keeps only several parameters that can generate the entire longitude and latitude values. This can save much space, but this also makes it difficult to retrieve those values unless HDF-EOS API is used.
The degree of difficulty of obtaining longitude and latitude values depends on the projection method.
The simplest projection is the geographic projection.
We will cover this case using a real AMSR-E AE_RnGd file from NSIDC.
You can download the file here.
To access a data field in a grid data, the enclosing HDF-EOS2 file should be opened using hdfgd, first.
Using the descriptor returned by that API, the grid should be opened.
In this example, we will access a grid named MonthlyRainTotal_GeoGrid in an HDF-EOS2 file named AMSR_E_L3_RainGrid_B05_200707.hdf.
hdfgd is the function which MATLAB provides to read data from an HDF-EOS2 grid file.
The first argument that it takes always reflects the task it is supposed to do.
In the figure above, open would open the file which is similar to GDopen API
from the HDF-EOS2 C library.
It will create a handle which is stored in file_id.
Using this file_id descriptor, the grid is opened in the next step using attach.
Then, data fields, dimensions and attributes in this grid can be accessed.
readfield.
The name of the data field being read is RrLandRain.
Note that grid_id, is the descriptor returned by attach.
rrland 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 rrland because MATLAB is a dynamically typed language.
After finishing reading all fields, one needs to close the grid using detach.
Similarly, if accessing the grid file is done, one needs to close the file using close.
grid_id is what attach returned, and file_id is what 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. For more information on how to retrieve latitude and longitude from an HDF-EOS2 grid, please check here .
The geographic projection, which is used by MonthlyRainTotal_GeoGrid, 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:
We need to retrieve those four types of information.
We will use gridinfo API to get the values for all the required variables.
xdimsize (72) and ydimsize (28) are the values for numX and numY, respectively. Also, upleft and lowright represents coordinates of upper left corner and lower right corner respectively in DMS format[3]. One needs to convert these coordinates in degree. After conversion, one can find these values as 0 degree, 70 degree, 360 degree and -70 degree. As the name of upleft and lowright imply, we can get 0 as leftX, 360 as rightX, 70 as upperY and -70 as lowerY.
Another factor that affects the longitude and latitude values is the Pixel Registration option.
Since this option can not be obtained.
HDFE_CENTER is assumed. HDFE_CENTER assigns 0.5 to
both offsetX and offsetY.
Now, we have enough information from the file.
Generating the correct plot with a geographic map still requires several adjustments for the arrangments of data and latitude and longitude.
We need to pass the latitude, longitude and rrland values to the contour function
in order to plot the data. To ensure the plot to be aligned with the world map generated by
the geoshow function, the first argument of contour must be longitude and the second
argument must be latitude. According to Figure 5,
the size of rrland is [72 28], which implies that the first dimension is xdim(lon)
and the second dimension is ydim(lat) in MATLAB.
This is also required by the contour fuction, which needs lon to be the first argument
and lat to be the second argument to ensure the plot to be overlaid with the world map correctly.
However, the fastest changing dimension of rrland is along the same longitude line. The size is 72.
This can be observed by Figure 7.
rrland is really stored as [28 72] in the file.
Therefore a transpose of the rrland is needed to generate the correct plot.
Displaying the world map with the geoshow function also requires the longitude range
starting from -180 degree and ending at 180 degree.
So the longitude needs to be translated from '2.5 degree to 357.5 degree' to '-177.5 degree to 177.5 degree'.
The latitude also needs to be translated from "north to south"(decreasing of the latitude)
to "south to north"(increasing of the latitude) with the geoshow function.
Accordingly, the latitude needs to be translated from '67.5 degree to -67.5 degree' to '-67.5 degree to 67.5 degree'.
Also the lower left corner in MATLAB is treated as the origin of the coordinate.
However, as we can see from the information obtained by gridinfo,
the origin is defined as the upper left corner in the file.
Hence, we need to flip rrland before passing it to the contour function.
The code section is listed below.
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 contour function, users may need to refer to
MATLAB's detailed document for more information about this function.
You can see the complete code from here. Figure 10 shows the result.
Unlike grid data, swath data contains geo-location fields such as Longitude and Latitude. This allows users to draw plots without the extra step to calculate longitude and latitude values; these values can be read from geo-location fields. However, swath data may contain dimension maps that allow the size of geo-location fields smaller or larger than that of the related data fields. If this is the case, retrieving longitude and latitude values becomes complicated. We will not cover the case when dimension maps are used to handle a swath.
We will use one swath file from AMSR-E/Aqua L2A Global Swath Spatially-Resampled Brightness Temperatures product. You can get one sample file here.
The step to open, read and close is similar.
hdfsw is used instead of hdfgd.
As the function name implies, functions that contain sw are for swath,
and the functions that contain gd are for grid.
The above code reads one data field 23.8H_Approx._Res.3_TB_(not-resampled)
and two geo-location fields, Latitude and Longitude.
Usually, a swath data contains both geo-location fields, and they provide longitude values and latitude values.
As the readfield function used to handle a grid file,
the readfield function for the swath file in the above piece of code stores
the data elements in the variable data.
Now that lat and lon contain the latitude and longitude values respectively,
one can pass those variables to contour function as the following arguments:
geoshow creates the world map,
with land area colored white.
The fifth and the sixth parameters of the geoshow function are FaceAlpha and 0.
These parameters imply that FaceAlpha value is set to zero.
This will clear the land areas of the map where the data overlaps with the map so that the plot of the data
will show up in those areas.
In this example, the first and second arguments of contour
are two-dimensional because both Longitude and Latitude are two-dimensional.
Then, each element of data is mapped to each element of lat and lon.
You can see the complete code here. Figure 13 shows the result.
