JWB
JWB

Reputation: 198

How to convert NASA GPM data from HDF5 to GeoTIFF using python or gdal?

I'm trying to convert an HDF5 file to GeoTIFF, but I am having trouble. It's specifically an HDF5 file from this source, I have had success using gdal_translate on other HDF5 files. An example file can be found here from NASA's ArthurHou Server (you might need to make a free account with NASA to download the data) or also here from a google drive. It's a pretty standard HDF5 file. It opens fine with Panoply but trying to convert to GeoTIFF or opening in QGIS is weird. It should look like a satellite swath but it looks like this

enter image description here

Here is my sample code using python and gdalwarp. I grab the GCPs to get the extent then use the extent in a gdalwarp command. I've also tried rasterio and am having no luck. I keep getting the same weird looking data. Any help would be appreciated.

subdataset_path = "//FS/CSF/flagShallowRain"

try:
    result = subprocess.run(
        ["gdalinfo", f"HDF5:{input_file}:{subdataset_path}"],
        capture_output=True,
        text=True,
        check=True,
    )
    output = result.stdout
except subprocess.CalledProcessError as e:
    print("Error running gdalinfo:", e)
    print("Output:", e.output)
    exit(1)

gcps = []
gcp_pattern = r'\((\d+\.\d+),(\d+\.\d+)\) -> \(([\d\.-]+),([\d\.-]+),\d+\)'

matches = re.findall(gcp_pattern, output)

gcps = []
for match in matches:
    pixel, line, lon, lat = map(float, match)
    gcps.append({"pixel": pixel, "line": line, "lon": lon, "lat": lat})

xmin = min(item['lon'] for item in gcps)
xmax = max(item['lon'] for item in gcps)
ymin = min(item['lat'] for item in gcps)
ymax = max(item['lat'] for item in gcps)

gdal_warp = f'gdalwarp -t_srs EPSG:4326 -te {xmin} {ymin} {xmax} {ymax} {input_filepath} {output_file}'

Upvotes: 2

Views: 145

Answers (0)

Related Questions