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

ET_Tif_Conversion

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.

• Vince: Sinusoidal projections can be one-way. Please edit the question to contain the exact commands you have tried.
• user3591836: I've added two things I tried, though they were mostly guessing.
• Jose: The answer you've got seems fairly complete to me, but you might want to read about how to go about doing this here
• user3591836: @AndreJ - I've had problems with that in other datasets. How have you overcome that?
• user3591836: This one, in particular: earlywarning.usgs.gov/fews/datadownloads/Global/ETa%20Anomaly Problems with Hawaii and New Zealand.

Answer 1

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):

enter image description here

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

• user3591836: Thanks a lot - I'll try this and comment again. But the MOD11 product is in sinusoidal projection: lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/… (in reply to your comment above).
• AndreJ: MOD11 is, but not the earlywarning you linked above. Maybe you inserted the wrong link.
• user3591836: yes sorry - that's true. I just meant I have problems with the 180° meridian in that dataset. When I convert the Tiff file to MBTiles, New Zealand disappears in the 2 most zoomed out levels.
• AndreJ: I think that is a different subject to ask for.
• user3591836: Can I just ask why you chose the value 0.01 for the -tr parameter? My colleague thinks .00836 is the right value to keep the pixels at 1000 meters.

Answer 2

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 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

• BenP: I have deleted from the GE question. The fact of the matter is that this answer is the simplest solution, using the most common tool GDAL. Numerous students of mine were unable to find a suitable answer hence the cross posting to related questions. I feel having a answer that links to a solution using the academic and industrial standard tooling is acceptable within the spirit of SE.
• senek: It worked for me. A great way to convert MODIS hdf files to tiff images on Python