MODIS REST API transform data to GeoTIFF with Python

Question

I'm trying to transform the response from a request with land surface temperature data from MODIS API to a GeoTIFF:

https://modis.ornl.gov/rst/api/v1/MOD11A2/subset?latitude=0&longitude=0&band=LST_Day_1km&startDate=A2001001&endDate=A2001001&kmAboveBelow=1&kmLeftRight=1

This request returns the following json:

{
 "xllcorner": "-9833349.10", 
 "yllcorner": "1575263.27", 
 "cellsize": 926.6254330558338, 
 "nrows": 3, 
 "ncols": 3, 
 "band": "LST_Day_1km", 
 "units": "Kelvin", 
 "scale": "0.02", 
 "latitude": 14.182305668542996, 
 "longitude": -91.20138157154014, 
 "subset": [
 {
 "modis_date": "A2020001", 
 "calendar_date": "2020-01-01", 
 "band": "LST_Day_1km", 
 "tile": "h09v07", 
 "proc_date": "2020010222551", 
 "data": [15263, 15232, 15219, 15234, 15217, 15219, 15224, 15227, 15299]
 }
 ]
}

I have a function in Python that receives the JSON from before as input and creates a GeoTIFF:

from osgeo import gdal, osr
def createGeoTif(data):
 values = np.array(data['subset'][0]['data'])
 values.shape=(data['nrows'],data['ncols'])
 values = (values * float(data['scale']))-273.15
 drv = gdal.GetDriverByName("GTiff")
 ds = drv.Create('./test.tif',data['ncols'],data['nrows'],1,gdal.GDT_Float32)
 xllcorner = float(data['xllcorner'])
 yllcorner = float(data['yllcorner'])
 x_resolution=y_resolution = data['cellsize']
 ds.SetGeoTransform([xllcorner, x_resolution, 0, yllcorner, 0, y_resolution])
 epsg = 4326
 srs = osr.SpatialReference()
 srs.ImportFromEPSG(epsg)
 dest_wkt = srs.ExportToWkt()
 ds.SetProjection(dest_wkt)
 ds.GetRasterBand(1).WriteArray(values)

The problem is that the GeoTIFF is placed far away where it should be. The problem I think is that the request returns the xllcorner and yllcorner in MODIS sinusoidal projection but I can't find in Google to which EPSG number it belongs so I just need to change this part of the code to fix it:

 epsg = 4326 #Modis sinusoidal code?
 srs = osr.SpatialReference()
 srs.ImportFromEPSG(epsg)
 dest_wkt = srs.ExportToWkt()
 ds.SetProjection(dest_wkt)

From modis api doc:

• xllcorner --- MODIS sinusoidal x-coordinate of the lower left corner of the subset
• yllcorner --- MODIS sinusoidal y-coordinate of the lower left corner of the subset

Answer 1

EPSG 4326 is lat lon coordinates.

The EPSG you should use for Modis sinusoidal projection is this one

Example code:

proj4 = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs "
srs = osr.SpatialReference()
srs.ImportFromProj4(proj4)
dest_wkt = srs.ExportToWkt()
ds.SetProjection(dest_wkt)