Tuesday, 4 December 2018

Reprojecting MODIS-Sphere ESRI GRID file to WGS84 using QGIS/GDAL


I have an ESRI GRID file projected in MODIS-Sphere (here's the dataset http://urban.yale.edu/sites/default/files/data/setoguneralphutyra_uf.zip) that I need to convert to a WGS84 GeoTIFF.


I tried Raster > Projections > Warp (Reproject) and got this nonsense: QGIS_Output



I then thought I'd use the MODIS reprojection tool or HEG to convert it, but they only accept HDF files and using QGIS' warp or translate functions to convert this file to HDF throwns an error. This guy was trying to do the same thing in Arc, no joy: MODIS-Sphere to WGS84 coordinate system in ArcGIS


I looked for clues in the metadata, and the file seems a little goofy; gdalinfo yields:


Driver: AIG/Arc/Info Binary Grid
Files: seto_uf
seto_uf.aux
seto_uf/.aux.xml
seto_uf/.DS_Store
seto_uf/dblbnd.adf
seto_uf/hdr.adf
seto_uf/log

seto_uf/metadata.xml
seto_uf/sta.adf
seto_uf/vat.adf
seto_uf/w001001.adf
seto_uf/w001001.tif
seto_uf/w001001.tif.aux.xml
seto_uf/w001001x.adf
Size is 7468, 2911
Coordinate System is `'
Origin = (-17345341.010356999933720,8423784.764957003295422)

Pixel Size = (5000.000000000000000,-5000.000000000000000)
Corner Coordinates:
Upper Left (-17345341.010, 8423784.765)
Lower Left (-17345341.010,-6131215.235)
Upper Right (19994658.990, 8423784.765)
Lower Right (19994658.990,-6131215.235)
Center ( 1324658.990, 1146284.765)
Band 1 Block=256x4 Type=Byte, ColorInterp=Undefined
Description = Seto_UF
Min=0.000 Max=101.000

NoData Value=255

So it can't determine the coordinate system.


There's a .tif included in the directory, but gdalinfo on that shows:


Driver: GTiff/GeoTIFF
Files: w001001.tif
w001001.tif.aux.xml
Size is 7468, 2911
Coordinate System is:
GEOGCS["Unspecified datum based upon the GRS 1980 Authalic Sphere",

DATUM["unknown",
SPHEROID["unnamed",6371007,1]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4047"]]
Origin = (-17345341.010356999933720,8423784.764957003295422)
Pixel Size = (5000.000000000000000,-5000.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:

INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-17345341.010, 8423784.765)
Lower Left (-17345341.010,-6131215.235)
Upper Right (19994658.990, 8423784.765)
Lower Right (19994658.990,-6131215.235)
Center ( 1324658.990, 1146284.765)

So it says the units are degrees but those huge corner coordinates and origin values don't agree.


Any ideas?



My setup:
OSX 10.9.4
QGIS 2.6.1 Brighton
GDAL 1.11.2 (apparently there are newer versions but this is the latest you can get on OSX without having to compile it yourself)



Answer



The Readme.docx inside the zip file gives some hints on the projection, but maybe misleading.


The raster is in an interrupted Goode Homolosine projection, based on the MODIS sphere.


So I created the following custom CRS:


+proj=igh +R=6371007 +no_defs


reprojected a Natural Earth shapefile into it, and it fits to the raster data with the same CRS assigned (not reprojected!):


enter image description here


So once you assign the custom CRS to the raster, you can set the project CRS on-the-fly to Web Mercator or something else:


enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...