0

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_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

0

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.