I'm trying to transform the response from a request with land surface temperature data from MODIS API to a GeoTIFF:
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