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
I haven't really gotten anywhere but this is the code I have so far
from osgeo import gdal
file_name='MCD19A2.A2010001.h11v11.006.2018047060424.hdf’
file_template = 'HDF4_EOS:EOS_GRID:"%s":grid1km:%s'
geojson=“polygon.geojson”
## 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
e=g_aod_warped.GetMetadata()['EASTBOUNDINGCOORDINATE']
w=g_aod_warped.GetMetadata()['WESTBOUNDINGCOORDINATE']
n=g_aod_warped.GetMetadata()['NORTHBOUNDINGCOORDINATE']
s=g_aod_warped.GetMetadata()['SOUTHBOUNDINGCOORDINATE']
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
{"type":"FeatureCollection","features":[{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[[[-70.3,-24.5],[-69.45,-24.5],[-69.45,-23.7],[-70.3,-23.7],[-70.3,-24.5]]]}}]}
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
