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.

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

3 Answers 3


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


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.

  • 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). May 31, 2016 at 4:06
  • MOD11 is, but not the earlywarning you linked above. Maybe you inserted the wrong link.
    – AndreJ
    May 31, 2016 at 5:46
  • 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. May 31, 2016 at 10:23
  • I think that is a different subject to ask for.
    – AndreJ
    May 31, 2016 at 14:52
  • 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. Jun 22, 2016 at 21:48

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
  • 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.
    – BenP
    May 24, 2021 at 16:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.