Adam,

It looks like you will not be able to use my HDF-EOS to Points tool to convert these files. The details are a little complicated.

Let’s start with the full-size file. As you mentioned, this file is very large. Each 2D variable (e.g. Rrs_413) has 12161 x 4481 cells. The tool expects to find latitude and longitude variables that are the same size, 12161 x 4481, giving the center coordinates of each cell. It uses these to create the points. But the latitude and longitude variables are only 12161 x 561, so the tool fails.

So what’s up with that? The problem is that HDF-EOS has several different ways for georeferencing grids. One scheme, designed to save space in the file, allows coordinates to be given for only every few cells, and requires the reader of the data to implement an algorithm to compute the values of the coordinates in between. So, rather than giving coordinates for all 4481 cells in that dimension, coordinates are only given for 561 cells, and the reader must determine the values for the rest.

This scheme is very complicated and I did not provide support for it in my tool, which I wrote in an afternoon to help someone out. Thus, my tool will not work on files that use this georeferencing scheme. (Really, the HDF-EOS designers say you are not even supposed to read the HDF directly, as my tool is doing. Instead, programs are supposed to integrate the HDF-EOS programming library with themselves and call functions in that library to read the data. That is the long term solution, but also complicated from a programming point of view, compared to just reading the HDF directly, and not something I could implement in that afternoon.)

(If I may get on my soap box for a moment, I believe it was not a good decision to use the complicated georeferencing scheme for the data we are dealing with here, because the relative space that is saved is negligible. This is because there are a number of variables in the file that have dimensions 12161 x 4481. Shrinking latitude and longitude to 12161 x 561 will save a little space, but all those others will still be there at the full 12161 x 4481 resolution. The end size of the file might shrink a few percent, but that is it. Why make reading the file much more difficult just to save a few percent of space?)

Now, what about the smaller file? In that file, each 2D variable only has dimensions 2905 x 4481. And it turns out that the latitude and longitude variables DO have the same dimensions as the other 2D variables (2905 x 4481). Thus it looked to me like this file would be convertible by my tool. But for some reason, the latitudes and longitudes appear to be messed up, with extreme ranges. This can be seen when you extract the HDF header, in the File Attributes section:

          Attr48: Name = Northernmost Latitude

                   Type = 32-bit floating point

                    Count= 1

                   Value = 2063693824.000000

           Attr49: Name = Southernmost Latitude

                   Type = 32-bit floating point

                    Count= 1

                   Value = -2777125632.000000

           Attr50: Name = Easternmost Longitude

                   Type = 32-bit floating point

                    Count= 1

                   Value = 3412.605957

           Attr51: Name = Westernmost Longitude

                   Type = 32-bit floating point

                    Count= 1

                   Value = -179.999405

When I debugged my script, I observed those same extreme values. I am not sure what is wrong with the coordinates. My guess is that whatever subsetting tool you used either stored bogus coordinates or it is using another complicated HDF-EOS georeferencing scheme that I do not understand and that my simple tool is not compatible with.

I guess at this point I can only recommend MRT or SeaDAS, unless you can figure out how to extract a subsetted file that has fully-populated, correct latitude and longitude values.

Best,

Jason

Adam,

An update:

The good news: the conversion only took 20 minutes on my server, instead of the 75 minutes I was expecting, because only 3.7 million points were inserted out of a possible 13 million. This presumably happened because most of the cells were No Data, in which case a point is not created in the output feature class.

The bad news: it looks like the x,y coordinates are screwed up, meaning there is still a bug in the script. I’m going to look at that now; it should only take a few minutes to fix.

Jason

Adam,

Thanks for the example file. After trying it out, I discovered that I accidentally sent you an earlier version of the tool. Here is a newer version. I tested it on your file and it appears to be working.

You recall me mentioning that it is pretty slow. Here are some measurements. The file contains about 13 million points. On my 18-month-old laptop, decently powerful for that time, with Arc 10.0 SP4, it is projected to take 6 hours. The bottleneck is that the Python script must go through ESRI’s geoprocessor API to insert the points into output feature class. ESRI improved this greatly in 10.1. For comparison: on a 4-month-old powerful server with an SSD hard disk and Arc 10.1, it is projected to take 75 minutes.

Obviously, the more you can trim down the area of interest, the fewer points you’ll have and the faster it will go.

If the slowness is a problem (and I certainly can understand if it is!) then you should try another tool, such as MRT. Given more time, I could write a better tool myself—it would use Python or MATLAB code to interpolate a raster directly from the non-uniformly-spaced array in the HDF, and go MUCH faster—but I am on paternity leave and already have some other problems I need to finish fixing right now. I hope you are able to make do with either the current tool or one of the other solutions.

CC: mget-help so anyone following gets the latest version of the tool.

Best,

Jason

Adam,

Sorry about that; it looks like a bug in my tool. Can you please email me the file M2010121163948.L2_FRS.bz2 so I can repro and fix it? I have not requested access to the MERIS data from NASA OceanColor so I cannot currently download it myself. Please be sure to send the bz2 file (so it is small enough to go through email) and please reply just to me, not mget-help, so we don’t spam everyone with the large file.

Thanks,

Jason

Hi Jason,

When running your hdf to point tool I get the error below.  Is there a method either in the script or tool that I need to change for the script to read the image values as integer instead of float or double? Or is that even the problem?

C:\seadas_shared\test\M2010121163948.L2_FRS.x.hdf\M2010121163948.L2_FRS.x.hdf.

TypeError: Cannot set the value of the Value field of this point of ArcGIS ShapeFile "C:\seadas_shared\test\Rrs_413.shp" to -32767.0 because the data type of that field is int16. To set int16 fields, you must provide an instance of <type 'int'> or <type 'long'>.

Finished inserting points: 0:00:00 elapsed, 0 points inserted, 0:00:00 per point.

Completed script ConvertHDFEOS2DSwathVariableToPoints...

Failed to execute (ConvertHDFEOS2DSwathVariableToPoints).

Adam,

Yeah, it’s too bad about Arc. I should say though that there might be some other way to do it, maybe using some Arc tools that were designed for LIDAR, which is a similar problem (taking points with arbitrary coordinates and conveniently interpolating a grid).

BTW, you do not need Linux to run GDAL. There are Windows builds available. You may be able to use GDAL to do the L2 HDF-EOS à raster interpolation. I have seen discussion of that on the gdal-dev list before.

The script is attached. Let me know if you have any problems. There is a toolbox for running it from ArcGIS.

I’m pretty sure MERIS uses HDF-EOS format so that MRT can probably handle it.

Best,

Jason

Hi Jason,

Thanks for the insight.  I had an inkling the cell size differences may be an issue but didn’t know Arc is unable to make the change. Goes to show Arc still isn’t the end all be all but integration with python helps.  I just installed a virtual linux OS on my computer, because I’d like to take advantage of batch processing of larger datasets in GDAL (still installing). I have SeaDas working fine just getting past the learning curve to figure out how the interpolation works will take a little time.  I would appreciate it if you could send me your Arc point to raster script to see how I can implement them both. So are you saying as well that the MRT can accept MERIS HDF as well? That’s probably the lightest way to do it.  I really appreciate your help

Best,

Adam Bausch

Adam,

Unfortunately I do not know of an easy way to work with L2 data in the current versions of ArcGIS (10.1 and earlier). Although L2 data typically come as 2-dimensional arrays in HDF or netCDF format, ArcGIS cannot currently understand how to georeference those arrays (unlike SeaDAS, which does understand them). The problem is that ArcGIS requires that all cells in an image have the same cell size. Unfortunately L2 data do not meet that requirement. L2 images are swaths or strips observed by the satellite. Cells near the center of the swath are spaced closer together than cells at the edge of the swath. ArcGIS does not know how to handle that.

When people encounter this problem, they naturally think it is an issue with projections—they suspect they are not using the right projection or the right parameters. Unfortunately that is not the problem. There are no projections you can use in ArcGIS that will allow you to georeference swath data correctly (at least that I know of).

To work around this you can either:

1.    Use SeaDAS to interpolate a raster image that does have a constant cell size, using a format that ArcGIS can read. GeoTIFF is a common choice here.

2.    If you are using MODIS data or other data in HDF-EOS format (I think MERIS will qualify), the MODIS Reprojection Tool (MRT) can do the same thing. The advantage of this over SeaDAS is that you can get a Windows version of the MRT. (SeaDAS requires you to run Linux in a virtual machine, at least the last time I checked, which can be complicated if you are used to Windows, and requires a lot of resources.)

3.    Use a geoprocessing tool I wrote to convert HDF-EOS L2 data to a point feature class, then use ArcGIS to interpolate a raster from that data using whatever projection and cell size you desire. This might be preferable if you prefer to work within ArcGIS and not learn a new tool like SeaDAS or MRT. But the HDF à Points à Raster workflow can be slow if you have many files to process. If you want to try this approach, let me know and I’ll send you the tool. (Eventually I’ll put this tool into MGET, once I can write a faster version of it.)

Hope that helps,

Jason

Hello,

I am a former Duke MEM and I’m trying to process some MERIS L2 FRS data obtained from the Ocean Color portal.  Before retrieving the data I had the imagery subset to the Gulf of Mexico region.  When viewing the image in Seadas (below)

the image footprint appears correct but when I use your MGET import HDF to ArcGIS raster tool, the image is shifted east and south when reprojected to NA Albers Equal Area Conic.  I might have made the erroneous assumption that the original Ocean Color data were in WGS1984. One question, is the input GCS correct? I filled in the script output parameters to those of the seadas band file (below). Assuming that the input cell size in decimal degrees is 0.002727 and the final output should be 300m.  Where am I going wrong? Also, I thought the L2 Rrs was a unitless value from 0-1 but from the band info that isn’t the case.  I appreciate any assistance you can lend to set me straight.

All the Best,

Adam Bausch

Adam Bausch
SENIOR STAFF SCIENTIST
CARDNO ENTRIX



Phone (+1) 813-664-4500  Fax (+1) 813-664-0440  Direct (+1) 813-367-0977  Mobile (+1) 919-428-0343  
Address 3905 Crescent Park Drive, Riverview, FL 33578 USA
Email Web www.cardno.com - www.cardnoentrix.com

This email and its attachments may contain confidential and/or privileged information for the sole use of the intended recipient(s). All electronically supplied data must be checked against an applicable hardcopy version which shall be the only document which Cardno warrants accuracy. If you are not the intended recipient, any use, distribution or copying of the information contained in this email and its attachments is strictly prohibited. If you have received this email in error, please email the sender by replying to this message and immediately delete and destroy any copies of this email and any attachments. The views or opinions expressed are the author's own and may not reflect the views or opinions of Cardno.