Earthdata Forum

Dear Support,
I'm a user of VIIRS Snow Cover products. I download L3 products (VNP10A1) with earthaccess. These data are distributed in h5 format and are not "natively" georeferenced, at least for common GIS softwares like GDAL or QGIS. In the user guide that I found on NSDIC website - VIIRS/[NPP|JPSS1] Snow Cover Daily L3 Global 375m SIN Grid, Version 2 taht I found here https://nsidc.org/data/vnp10a1/versions/2 - there is the information regarding the geometry of the product, so that an user can apply georeferencing from the information contained in the products (XDim, YDim coordinates). In this document the CRS information is given with a PROJ4 string. However, in turned out in my experience that this string is wrong and the MODIS CRS PROJ4 string should be used to have correct georeferencing.

User guide PROJ4 string:
+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
MODIS PROJ4 string:
+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs

Is it normal that VIIRS snow cover product georeferencing cannot be read by GDAL/QGIS ? What about the CRS, is my finding correct? (can be helpful for other users if corrected in new versions of the documentation)
Thank You for your amazing work.

Nicola Imperatore

Dear Nicola Imperatore,

Thank you for your question and observations regarding VIIRS Snow Cover L3 CRS/PROJ4 string information found in the User Guide for these products. The NSIDC is taking a look into your inquiry and will follow up shortly with more information.

Kind Regards,
Gail

Hi Nicola

You are correct, the coordinate reference system (CRS) for the VIIRS snow cover data is `+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs`

We will be changing the user guides to reflect this.

You can load data from HDF5 files into QGIS but you will have to define the CRS and grid definition manually. GDAL is a useful tool for this. MyHDF5 is also a useful tool https://myhdf5.hdfgroup.org/help.

Using MyHDF5 you can see a number of groups. Navigate to the /HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data Fields. You will see a Projection data variable.

Using GDAL, you can view the projection information
$ gdalinfo HDF5:"data/VNP10A1.A2024168.h15v02.001.2024169145053.h5"://HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/Projection

This returns a listing, including the CRS in Well Known Text (WKT) format.

To access a data variable in QGIS, you need to find the path to the data variable. You can use gdalinfo for this.

$ gdalinfo data/VNP10A1.A2024168.h15v02.001.2024169145053.h5

This produces a long listing. The relevant section is
<snip>
Subdatasets:
SUBDATASET_1_NAME=HDF5:"data/VNP10A1.A2024168.h15v02.001.2024169145053.h5"://HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/Algorithm_bit_flags_QA
SUBDATASET_1_DESC=[3000x3000] //HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/Algorithm_bit_flags_QA (8-bit unsigned character)
SUBDATASET_2_NAME=HDF5:"data/VNP10A1.A2024168.h15v02.001.2024169145053.h5"://HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/Basic_QA
SUBDATASET_2_DESC=[3000x3000] //HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/Basic_QA (8-bit unsigned character)
SUBDATASET_3_NAME=HDF5:"data/VNP10A1.A2024168.h15v02.001.2024169145053.h5"://HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/NDSI
SUBDATASET_3_DESC=[3000x3000] //HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/NDSI (16-bit integer)
SUBDATASET_4_NAME=HDF5:"data/VNP10A1.A2024168.h15v02.001.2024169145053.h5"://HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/NDSI_Snow_Cover
SUBDATASET_4_DESC=[3000x3000] //HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/NDSI_Snow_Cover (8-bit unsigned character)
SUBDATASET_5_NAME=HDF5:"data/VNP10A1.A2024168.h15v02.001.2024169145053.h5"://HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/granule_pnt
SUBDATASET_5_DESC=[3000x3000] //HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/granule_pnt (8-bit unsigned character)

<snip>

I'll assume you are interested in the NDSI_Snow_Cover variable.

In QGIS open the add_layer dialog and copy and paste the full SUBDATASET name for NDSI_Snow_Cover:

HDF5:"data/VNP10A1.A2024168.h15v02.001.2024169145053.h5"://HDFEOS/GRIDS/NPP_Grid_IMG_2D/Data_Fields/NDSI_Snow_Cover

and then add the layer. This layer will not have projection information. You will need to add this manually.

See here https://docs.qgis.org/3.40/en/docs/user_manual/working_with_projections/working_with_projections.html#custom-coordinate-reference-system

Go to Settings > Custom Projections to create a custom projection. Give the new projection a name, select Format as WKT and copy and paste the WKT for the Sinusoidal projection into the parameters box.

Now click on the projection icon in the lower left corner of the map and select and apply the projection you just defined.

You can set the extent of the data in the Layer Properties > Metadata > Extent. You will need to get the x and y extents from the XDim and YDim variables which you can view in MyHDF5

I hope that helps.

Andy Barrett
NSIDC

Hi Nicola,

All that Andy wrote is correct! As an alternative to the instructions Andy provided, with more recent versions of QGIS, it is possible to use the GUI if you prefer and drag the .h5 file into the Layers panel on the left. You'll then see the options for the parameters in the data that you can load (see attached figure with the pink arrow suggesting the drag done).

It will load with an "invalid projection" and you'll have to assign the custom CRS as Andy described, but the extent information will be available in the Properties > Information menu off the layer in the layers panel (right click on the layer to access this menu). This alleviates the need for MyHDF5.

Cheers,
Diane

Dear Gail, Andy and Diane,
Thank you for your readiness and your complete answers. The tipps you gave me are very useful.
In my case, I had a long time series to process and at the same time I wanted to be able to inspect the data puncutally on QGIS for effective comparison with other satellite data. I ended up adding automazing the georeferencing process in Python. I leave the code snippet here in case other users are interested.

import xarray as xr
import rioxarray

# Read Grid information using xarray
grid_data = xr.open_dataset(viirs_l3_filepath, group="HDFEOS/GRIDS/VIIRS_Grid_IMG_2D", engine="netcdf4")
bin_size = xr.open_dataset(viirs_l3_filepath, engine="netcdf4").attrs["CharacteristicBinSize"]

# Compute x and y coordinates
# ! The grid information in VIIRS snoow cover product refers to upper left corner of the pixel
x0 = grid_data.coords["XDim"][0]
y0 = grid_data.coords["YDim"][0]
width=len(grid_data.coords["XDim"])
height=len(grid_data.coords["YDim"])
xmin = x0 + bin_size / 2
xmax = xmin + (width - 1) * bin_size
ymin = ymax - (height - 1) * bin_size
ymax = y0 - bin_size/2
xcoords = np.linspace(xmin, xmax, width)
ycoords = np.linspace(ymax, ymin, height)

# Read snow cover data
ndsi_snow_cover = xr.open_dataset(viirs_l3_filepath, group="HDFEOS/GRIDS/VIIRS_Grid_IMG_2D/Data Fields", engine="netcdf4").data_vars["NDSI_Snow_Cover"]
# Apply coordinates
ndsi_snow_cover = ndsi_snow_cover.rename({"XDim": 'x', "YDim": 'y'}).assign_coords(coords={'y: ycoords, 'x': xcoords})

# Force CRS using rioxarray
modis_crs = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"
ndsi_snow_cover_georefrenced = ndsi_snow_cover.rio.write_crs(modis_crs).rio.write_coordinate_system()

This looks to be working in my case.
Best Regards,

Nicola

I've made a note of this and will try to get this fixed in a future GDAL, this is presumably a common way to georeference with the XDim/YDim coord values.

Just generally, I would suggest reporting examples like this directly as issues on the repo, or ask on gdal-dev, there's no end to the variations that can be found in array formats and GDAL is a well supported and responsive project. (I will do this at some point but will also do extra pre-work in case I can make the required changes in a PR myself).