2

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
2

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)

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.