5

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.


Edit: 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.


Edit: 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.


Comment: 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.

  • 1
    Sinusoidal projections can be one-way. Please edit the question to contain the exact commands you have tried. – Vince May 23 '16 at 10:57
  • I've added two things I tried, though they were mostly guessing. – user3591836 May 23 '16 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 '16 at 17:01
  • @AndreJ - I've had problems with that in other datasets. How have you overcome that? – user3591836 May 27 '16 at 1:25
  • This one, in particular: earlywarning.usgs.gov/fews/datadownloads/Global/ETa%20Anomaly Problems with Hawaii and New Zealand. – user3591836 May 27 '16 at 1:36
3

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

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.