I have a Python code that reads HDF4 files for different satellite products. I am having issues obtaining the underlying grid that it uses.

In panoply I can see the variables XDim and YDim that have the points in the sinusoidal grid but I cant seem to find a way to extract them with GDAL. I need to get the mesh in lat-lon for the whole tile but also for a smaller polygon

enter image description here

I haven't really gotten anywhere but this is the code I have so far

from osgeo import gdal

file_template = 'HDF4_EOS:EOS_GRID:"%s":grid1km:%s'

## Open HDF file
g = gdal.Open(file_name)

## Retrieve original data
this_file = file_template%(file_name, "Optical_Depth_055")
g_aod = gdal.Open(this_file)
proj_aod = g_aod.GetProjection() 

## Warp to geojson
g_aod_warped = gdal.Warp('', this_file, format = 'MEM', dstNodata = None, outputType = gdal.GDT_Float32, cutlineDSName = geojson, cropToCutline = True, resampleAlg = 0, dstSRS = 'EPSG:4326')
proj_aod_warped = g_aod.GetProjection() 

## Bounds

All the responses I've found on how to extract the grid assume an equal interpolation between the bounds but that wouldn't be true all the time.

Is there a way to get the grid with gdal?

I am seriously considering using pyhdf to retrieve that data, turn it into a geopandas and reproject it like that, that way I can also get the intersection with the polygon in a way that I know what's happening since If I look at the bounds variables they are outside the bounds of the geojson.

The geojson file is


More information:

>>> g_aod.GetProjection()                                                                                                                               
'PROJCS["unnamed", GEOGCS["Unknown datum based upon the custom spheroid", DATUM["Not specified (based on custom spheroid)", SPHEROID["Custom spheroid",637 1007.181,0]], PRIMEM["Greenwich",0],UNIT["degree”,  0.0174532925199433, AUTHORITY["EPSG", "9122"]]], PROJECTION["Sinusoidal"], PARAMETER["longitude_of_center" ,0], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["Meter",1], AXIS["Easting",EAST], AXIS["Northing",NORTH]]'   

>>> g_aod_warped.GetProjection()                                                                                                                                                  
'GEOGCS["WGS 84",DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AXIS["Latitude",NORTH], AXIS["Longitude",EAST], AUTHORITY["EPSG","4326"] ]'   

Python 3.10 GDAL 3.2.1


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.