Convert MODIS HDF file in sinusoidal projection into GeoTIFF using Python GDAL

Question

I'm looking for a way to convert MODIS LST data from HDF to Tiff format, and reproject from sinusoidal to WGS84 - all using Python GDAL.

An example file is here (I only care about band 0) : http://e4ftl01.cr.usgs.gov/MOLT/MOD11A1.005/2004.03.31/MOD11A1.A2004091.h34v10.005.2007261231833.hdf

After converting all the tiles for one date, I want to use gdalmerge to merge all them into a Tiff file for the entire globe.

I've tried various combinations of gdalwarp and gdal_translate, but I think I'm not specifying the parameters properly.

Here is something I tried.

gdal_translate -of GTiff {inputfile} {outputfile}

Together with one of the following (neither worked):

gdalwarp -of GTIFF -s_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' -r cubic -t_srs '+proj=longlat +datum=WGS84 +no_defs' {inputfile} {outputfile}


gdalwarp -s_srs '+proj=sinu +wktext' -t_srs EPSG:32719 -r cubic  {inputfile} {outputfile}

I'm new to GIS, though, and have been trying to figure it out from other stackexchange examples.

For anyone familiar with ArcGIS, here is a script that achieves what I'm hoping to do, but using ArcPy. It was written for MODIS ET data. Example here: ftp://ftp.ntsg.umt.edu/pub/MODIS/NTSG_Products/MOD16/MOD16A2.105_MERRAGMAO/Y2000/D193/MOD16A2.A2000193.h08v05.105.2013121202542.hdf

I want to replicate this using GDAL.

---

In the end I used the solution of Loïc Dutrieux for tiles away from the dateline, and for those near the dateline, I use the clipping solution of AndreJ.

Answer 1

You can either leave the -s_srs empty or use +proj=sinu +R=6371007.181 +nadgrids=@null +wktext

gdalwarp -of GTIFF -s_srs '+proj=sinu +R=6371007.181 +nadgrids=@null +wktext' -r cubic -t_srs '+proj=longlat +datum=WGS84 +no_defs' {inputfile} {outputfile}

Note that inputfile must be the full subdataset (SDS) name, which you can get from gdalinfo. For example, HDF4_EOS:EOS_GRID:"MOD16A2.A2000193.h13v12.105.2013121204930.hdf":MOD_Grid_MOD16A2:ET_1km

Answer 2

You might get unexpected results because your dataset crosses the 180° meridian. As a consequence, the tile is squeezed around the globe when reprojected to WGS84. To avoid that, you have to cut the raster data at the +/- 179.99° meridian. The following batch works with pure GDAL:

gdal_translate HDF4_EOS:EOS_GRID:"MOD11A1.A2004091.h34v10.005.2007261231833.hdf":MODIS_Grid_Daily_1km_LST:LST_Day_1km LST.tif
gdaltindex -tileindex location LST-extent.shp LST.tif
ogr2ogr -segmentize 1000 LST-segmented.shp LST-extent.shp
ogr2ogr -t_srs EPSG:4326 -wrapdateline LST-extent-wgs84.shp LST-segmented.shp
ogr2ogr -clipsrc "POLYGON ((0 89.99, 179.99 89.99, 179.99 -89.99, 0 -89.99, 0 89.99))" LST-clipped.shp LST-extent-wgs84.shp
gdalwarp -overwrite -t_srs EPSG:4326 -tr 0.01 0.01 LST.tif LST-wgs84.tif
gdalwarp -overwrite -cutline LST-clipped.shp -crop_to_cutline -tr 0.01 0.01 LST-wgs84.tif LST-wgs84_clipped.tif

In short, I extract the LST band from the HDF to a tif file, then create a polygon from the extent, densify the lines to 1km interval, reproject it to WGS84 with a cut at the dateline, then clip the result to the Eastern hemisphere.

Then reproject the raster to WGS84, and clip it to the polygon.

With this result (the clip polygon is dashed):

The only thing that gets lost is the scale factor of 0.02 in the metadata.

Answer 3

For Python users, I have written Geospatial Analysis: obtaining and pre-processing OpenSource satellite data and a Gist explaining how to perform this translation in GDAL and Rasterio.

from osgeo import gdal
  
in_file = f"./modis.hdf" # raw MODIS HDF in sinusoid projection
out_file = f"./modis.tif"
    
# open dataset
dataset = gdal.Open(in_file,gdal.GA_ReadOnly)
subdataset =  gdal.Open(dataset.GetSubDatasets()[0][0], gdal.GA_ReadOnly)

# gdalwarp
kwargs = {'format': 'GTiff', 'dstSRS': 'EPSG:4326'}
ds = gdal.Warp(destNameOrDestDS=out_file,srcDSOrSrcDSTab=subdataset, **kwargs)
del ds