Open Georreferenced MODIS L2 atmospheric data in python

Question

i'm working with a MOD04 (atmospheric product) but i don't know how to georrefence the image in Python. I'll show an example of a code opening the Optical Thickness product

#==============================================================================
# MOD04 AOT PRODUCT
#==============================================================================
file_name04 = ("C:\MODIS\MOD04_L2.A2017112.1450.006.2017113014708.hdf")
MOD04 = gdal.Open(file_name04)
subdataset04=MOD04.GetSubDatasets()
meta04=MOD04.GetMetadata()
MODAOD = gdal.Open(subdataset04[63][0]) ###AOT PRODUCT
AOT = MODAOD.ReadAsArray()
meta04aot=MODAOD.GetMetadata()
AOT= AOT.astype(np.float32)
AOT[AOT == -9999] = np.nan
mult=float(meta04aot['scale_factor'])
add =float(meta04aot['add_offset'])
AOT = mult*(AOT-add)

The GeoTransform showed in the image notes that is not georrefenced:

 MODAOD.GetGeoTransform()
 >> (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
I try to set the boundaries using MODAOD.SetGeoTransform() with the north and west coordinates.
 XSize = MODAOD.RasterXSize
 YSize = MODAOD.RasterYSize
XSize
>> 135
YSize
>> 203
MOD04lat = meta04aot['GRINGPOINTLATITUDE.1']
MOD04lon = meta04aot['GRINGPOINTLONGITUDE.1']
MOD04lat,MOD04lon=MOD04lat.split(","),MOD04lon.split(",")
MOD04lat
>>: 
['-4.89502619820019',
 ' -7.98504122087105',
 ' -26.0833829669723',
 ' -22.6082546570904']
MOD04lon
>>: 
['-76.310463096777',
 ' -55.2650241843125',
 ' -58.3815019710174',
 ' -81.2664506382633']
north=np.float(meta04['NORTHBOUNDINGCOORDINATE'])
west =np.float(meta04['WESTBOUNDINGCOORDINATE'])
 north
>>: -4.92184878405504
west
>>: -81.2587263920377
pixY_size = (abs(float(MOD04lat[0])-float(MOD04lat[2]))/YSize)
pixX_size = (abs(float(MOD04lon[3])-float(MOD04lon[2]))/XSize)
pixY_size
>> 0.1043761417180892
pixX_size
>> 0.1695181382758955
geo_transform=(float(MOD04lon[3]),pixX_size,0,float(MOD04lat[0]),0,
(pixY_size*-1))
MODOAOD.SetGeoTransform(geo_transform)
MODOAOD.SetProjection('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,\
AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],\
PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]')
output=("C:\MODIS\MOD04_AOT.tif")

when i export the image i use this code

def save_raster ( output_name, dataset, raster_data, driver ,NaN_Value,arr_projection=None):
 if arr_projection is None:
 arr_projection=[]
 else:
 if str(type(arr_projection)) == "<class 'osgeo.gdal.Dataset'>":
 srs=arr_projection.GetProjectionRef ()
 else:
 srs = arr_projection
 """
 A function to save a 1-band raster using GDAL to the file indicated
 by ``output_name``. It requires a GDAL-accesible dataset to collect 
 the projection and geotransform.
 """
 # Open the reference dataset
 g = ( dataset )
 # Get the Geotransform vector
 if raster_data is False:
 raster_data=g.ReadAsArray()
 if type(raster_data) == tuple:
 raster_data = np.array(raster_data[0])
 if str(type(g)) == "<class 'osgeo.gdal.Dataset'>":
 geo_transform = g.GetGeoTransform ()
 x_size = g.RasterXSize # Raster xsize
 y_size = g.RasterYSize # Raster ysize
 srs = g.GetProjectionRef () # Projection
 elif str(type(g)) == "<class 'affine.Affine'>":
 geo_transform = (g[2],g[0],g[1],g[5],g[3],g[4])
 RastArr = raster_data
 x_size = int(RastArr.shape[1])
 y_size = int(RastArr.shape[0])
 #PROCESS RASTERIO NUMPY
 else:
 geo_transform = (g[1][2],g[1][0],g[1][1],g[1][5],g[1][3],g[1][4])
 RastArr = np.array(g[0])
 x_size = int(RastArr.shape[2])
 y_size = int(RastArr.shape[1])
 if raster_data.ndim > 2:
 raster_data=raster_data[0]
 NaN_rast = NaN_Value
 # raster_data[raster_data == NaN_rast] = 'NaN'
 raster_data[raster_data == NaN_rast] = np.NaN
 # Need a driver object. By default, we use GeoTIFF
 driver = gdal.GetDriverByName ( driver )
 dataset_out = driver.Create ( output_name, x_size, y_size, 1, \
 gdal.GDT_Float32 )
 dataset_out.SetGeoTransform ( geo_transform )
 dataset_out.SetProjection ( srs )
 dataset_out.GetRasterBand ( 1 ).WriteArray ( \
 raster_data.astype(np.float32) )
save_raster(output,MODAOD,AOT,"GTiff",300000)

but when i open the image exported with python vs the image opened with qgis show that is the boundaries are not well georeferenced, because the MODIS image is not exactly rectangular

A picture taken from the hdf opening QGIS (well georefenced) enter image description here

and the image exported with the python code

enter image description here

anyone know a pythonic solution?

Answer 1

Found a solution using gdal warp plus a thin-plate spline interpolation..

MOD04 = gdal.Open(file_name04)
MAOT=gdal.Open(MOD04.GetSubDatasets()[63][0]) #AOT EXAMPLE
gdal.Warp('C/directory/imagery.tif',MAOT,tps=True)

• Comrade Che: Does GDAL support the MODIS HDF files georeferencing? Is that python code? If so how the second line works? This space between variables is syntax error: MAOT gdal.Open(...), no?
• Italo Moletto Lobos: Mr. Che, i've missed the equal in the code line, so i fixed the code line. MAOT is a object created by the command gdal, open, so i added the equal ("=") between MAOT and gdal.Open