# "Good" HDF

## Data source

[MOD10A ](https://nsidc.org/data/myd10a1/versions/6)/ [MYD10A](https://nsidc.org/data/myd10a1/versions/6) - daily snow cover, identified using Normalized Difference Snow Index (NDSI). Files could be downloaded directly from NSIDC, or from NASA Earthdata. 

{% hint style="info" %}
For this section h15v14 tile (covering South Georgia) was used.

Files were downloaded from Earthdata with `wget`. For this reason, the folders structure is the following: `MOD10A1/n5eil01u.ecs.nsidc.org/DP4/<sensor code>/<product code>/<YYYY.MM.DD>/<filename>`
{% endhint %}

## QGIS

### Add data to the map view

This window appears after adding an hdf file to the map: sub-dataset selection is needed. NDSI snow cover data is stored in `'NDSI_Snow_Cover MOD_Grid_Snow_500m (8-bit unsigned integer)'` sub-dataset.

<figure><img src="https://1684505540-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2Ffw43Ebzk64frU0F55GJ4%2Fuploads%2FBWlV0ZwKiIIXQNWLjbfw%2Fimage.png?alt=media&#x26;token=1af422f4-dd77-40f0-bb73-66e830356b91" alt=""><figcaption><p>QGIS dialog window for sub-dataset selection</p></figcaption></figure>

### Check the georeference

Data is georeferenced correctly, though QGIS doesn't know the MODIS data CRS, as it is not in the EPSG catalog.

<figure><img src="https://1684505540-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2Ffw43Ebzk64frU0F55GJ4%2Fuploads%2FlZUBFBv0buUabeT3wWrK%2Fimage.png?alt=media&#x26;token=cb28db7d-df01-4453-b0aa-7e005ad645ce" alt=""><figcaption><p>MYD10A hdf file added to a map in EPSG:3031</p></figcaption></figure>

### Export

Anyway, it is possible to work with this file, for example to reproject it and export to GeoTIFF format.

<figure><img src="https://1684505540-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2Ffw43Ebzk64frU0F55GJ4%2Fuploads%2FJ5FWzbOETeT9P5a9r5pC%2Fimage.png?alt=media&#x26;token=9ac30ab7-9a49-4e5c-a145-52deea91f7c9" alt=""><figcaption><p>Exporting MODIS HDF to GeoTIFF and clipping to vector data extent</p></figcaption></figure>

## GDAL

What if you need to process this way not only 1 file but at least 100? GDAL command line tool could be a better solution for this. GDAL installation instructions are on the [first ](https://guides.geospatial.bas.ac.uk/working-with-scientific-data-formats-in-gis/intro-notes#gdal-installation)page.

### gdalinfo: view metadata

Navigate to your working folder, e.g. `cd /mnt/c/my_dir/`. Run gdalinfo command on one of the HDF files:

{% code overflow="wrap" %}

```bash
gdalinfo MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf
```

{% endcode %}

The output contains all the metadata, but for further work, you need only a correct sub-dataset name. You will find it at the end of the output (`SUBDATASET_1_NAME`):

```bash
Driver: HDF4/Hierarchical Data Format Release 4
Files: MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf
Size is 512, 512
Metadata:
  ...
Subdatasets:
  SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:NDSI_Snow_Cover
  SUBDATASET_1_DESC=[2400x2400] NDSI_Snow_Cover MOD_Grid_Snow_500m (8-bit unsigned integer)
  SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:NDSI_Snow_Cover_Basic_QA
  SUBDATASET_2_DESC=[2400x2400] NDSI_Snow_Cover_Basic_QA MOD_Grid_Snow_500m (8-bit unsigned integer)
  SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:NDSI_Snow_Cover_Algorithm_Flags_QA
  SUBDATASET_3_DESC=[2400x2400] NDSI_Snow_Cover_Algorithm_Flags_QA MOD_Grid_Snow_500m (8-bit unsigned integer)  SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:NDSI
  SUBDATASET_4_DESC=[2400x2400] NDSI MOD_Grid_Snow_500m (16-bit integer)
  SUBDATASET_5_NAME=HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:Snow_Albedo_Daily_Tile
  SUBDATASET_5_DESC=[2400x2400] Snow_Albedo_Daily_Tile MOD_Grid_Snow_500m (8-bit unsigned integer)
  SUBDATASET_6_NAME=HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:orbit_pnt
  SUBDATASET_6_DESC=[2400x2400] orbit_pnt MOD_Grid_Snow_500m (8-bit integer)
  SUBDATASET_7_NAME=HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:granule_pnt
  SUBDATASET_7_DESC=[2400x2400] granule_pnt MOD_Grid_Snow_500m (8-bit unsigned integer)
```

### gdalwarp: clip, reproject and save to GeoTiff&#x20;

[gdalwarp ](https://gdal.org/programs/gdalwarp.html)command is used to reproject data with the option to clip the image to a certain extent. The full expression is the following

{% code overflow="wrap" %}

```bash
gdalwarp HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:NDSI_Snow_Cover modis_clip_gdal.tif -t_srs epsg:3031 -cutline SouthGeorgiaExtent.shp -crop_to_cutline -tr 500 500
```

{% endcode %}

* The first parameter is an input file, the sub-dataset specified as it appears in gdalinfo output: `HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:NDSI_Snow_Cover`
* The second parameter is the output file, which will be created in the working directory: `modis_clip_gdal.tif`
* `-t_srs` tag sets output spatial reference (epsg code in this case)
* `-cutline` and `-crop_to_cutline` tags let cropping the image to the vector file extent
* `-tr` tag defines target raster resolution (x and y) in CRS specified by `-t_srs`

To be sure, that this file is processed correctly, add it to the map view:

<figure><img src="https://1684505540-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2Ffw43Ebzk64frU0F55GJ4%2Fuploads%2FMMByWq9HBbFIiWNpsf2l%2Fimage.png?alt=media&#x26;token=a84b3a75-a137-4374-8c9c-2bbde9ad2089" alt=""><figcaption><p>Processed file in the map view</p></figcaption></figure>

### Batch  gdalwarp

To apply the same command to all the downloaded files, we need to amend the previous expression:&#x20;

{% code overflow="wrap" %}

```bash
for i in MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOS*/M*D10A1.061/*/*.hdf
do
gdalwarp -t_srs epsg:3031 -cutline SouthGeorgiaExtent.shp -crop_to_cutline -te_srs epsg:3031 -tr 500 500 HDF4_EOS:EOS_GRID:"$i":MOD_Grid_Snow_500m:NDSI_Snow_Cover $i.tif
done
```

{% endcode %}

This command creates GeoTIFF files with the same names and the same directories as HDF files. To move them to another folder:

<pre class="language-bash"><code class="lang-bash"><strong>mkdir MOD10A1_tiff
</strong><strong>mv MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOS*/M*D10A1.061/*/*.tif MOD10A1_tiff/
</strong></code></pre>

## Metadata

You probably have already noticed, that the NDSI snow cover layer has 8-bit values from 0 to 250, but it is not clear what these values stand for. The answer could be found with QGIS or GDAL.

Open HDF layer properties in QGIS and in "More information" section find '`Key= 0-100=NDSI snow, 200=missing data, 201=no decision, 211=night, 237=inland water, 239=ocean, 250=cloud, 254=detector saturated, 255=fill`'.

The same output you can get if run `gdalinfo` on a specified sub-dataset:

{% code overflow="wrap" %}

```bash
gdalinfo HDF4_EOS:EOS_GRID:"MOD10A1/n5eil01u.ecs.nsidc.org/DP4/MOSA/MYD10A1.061/2023.01.07/MYD10A1.A2023007.h15v14.061.2023009201030.hdf":MOD_Grid_Snow_500m:NDSI_Snow_Cover
```

{% endcode %}