4

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?

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)
  • 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? – Mr. Che Jul 27 '18 at 6:51
  • 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 – Italo Moletto Lobos Jul 28 '18 at 20:03

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.