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)