net2cog should support HDF4 format

[jackiryan]

There are many legacy HDF4 format datasets on the Earthdata cloud, especially for MODIS. Here's an example: https://www.earthdata.nasa.gov/data/catalog/lpcloud-myd09gq-061
Currently, net2cog does not support this format and there is no alternative harmony service for reformatting these datasets. In the code, you can see that these files are skipped right now:
https://github.com/podaac/net2cog/blob/develop/net2cog/netcdf_convert.py#L453-L483

    if netcdf_file.endswith(('.nc', '.nc4', 'h5')):
        ...
    else:
        logger.info("Not a NetCDF file; Skipped file: %s", netcdf_file)
        output_files = []

[tloubrieu-jpl]

Is HDF5 already supported ?

[jackiryan]

It should be, since HDF5 is effectively the same as netcdf 4, which is tested an supported. But we do not have specific tests for h5 at this time.

[D-Auty]

Cloud OPeNDAP is a potential prefix service in the case of HDF4 inputs, though I'm not sure OPeNDAP will generate sufficiently robust NetCDF outputs for Net2Cog to work correctly.

xarray itself is unlikely to be able to handle HDF4, but rasterio may be able to pre-process the hdf4 input file.

Rasterio is not necessarily ideal in the general hdf4-to-NetCDF conversion, but it might work remarkably well for eventual conversion to GeoTIFF. Rasterio's flaws in the general case - flattening of the variable hierarchy, mandatory inclusion of an empty "Band" dimension, etc. I suspect do matter when converting to GeoTIFF.

[jackiryan]

@D-Auty We'd have to add in the pyhdf module as a dependency and detect when the input to net2cog is hdf. I think this is a better solution since it consolidates reformatting functionality into the net2cog service. That said, it doesn't fix any challenges that may arise further up a service chain. For example, if we had an L2 HDF4 granule that needed swath projection and reformatting to PNG, the swath projection service also needs to be made HDF4 compatible.

Having a dedicated HDF4->nc microservice is flexible, but creates potential data duplication issues where we're storing more intermediate files. This will come up most often with MODIS collections, so the data volume is quite high (100s of GB/day)

[D-Auty]

This seems like a simple update, but ultimately has proved disappointing. The end-result is available as a GeoTIFF file, but it is not properly geolocated data.

I am surprised, as the xarray.to_netcdf output is geolocated. The .nc file is not great as a capture of the HDF data, from my perspective, but it is properly geolocated. Yet when converted from there to GeoTIFF (COG), something is lost.

This might be a rasterio, rioxarray or rio_cogeo issue, though, and we might be able to push through a fix there, which would make this approach viable again.

Note, re. LP-DAAC’s use-cases - I suspect this still will not work for ASTER, a collection which is a current topic for ECS-Decommissioning and potential candidate for support in the cloud. The ASTER data fails for other reasons having to do with the quirkiness of the ASTER data in HDF (HDFEOS2) format.

I added the following to NetCDF, netcdf_converter function (~ line 465)

elif netcdf_file.endswith(('.hdf', '.he4', 'h4')):
    logger.info("Reading %s", basename(netcdf_file))

    input_datatree = xr.DataTree(
        rioxarray.open_rasterio(
            netcdf_file,
            # masked=True,
            # mask_and_scale=True,
            mode='r',
    ))
    input_datatree.to_netcdf(output_directory+'/tmp.nc')
        # just for reference - the input file written as NetCDF data
else:
    logger.info("Not a NetCDF file; Skipped file: %s", netcdf_file)
    output_files = []

if input_datatree:

# ... continue with processing var_list

This shows how easy the update seems to be.

N.B. This might not be relevant, given that the end-result is still not geolocated - but I had this issue in other efforts to create GeoTIFF from HDF files. Rio_xarray introduces a single element “band” dimension that is not present in the native HDF file. This complicates the GeoTIFF generation. I only succeeded in these other cases when I removed that band dimension.

In _write_cogtiff method (~ line 110)

variable = nc_xarray[variable_path]
for dim in variable.dims:
    if dim == "band" and variable.sizes[dim] == 1:
        variable.drop_vars(dim)
nc_xarray[variable_path].rio.to_raster(temp_file_name)