Search code examples
pythongdalprojectiongeotiffhdf

Converting HDF files to GeoTiff files using Python's GDAL

I want, convert the file named "MAM35S0.A2008010.0000.002.2008010173218.hdf"(MODIS Aqua product) to GTiff format.

The datasheet is "HDF4_EOS:EOS_SWATH:"MAM35S0.A2008010.0000.002.2008010173218.hdf":mod35:Cloud_Mask", but normal conversion was not possible using the API, so I used the command.

from osgeo      import gdal

import numpy                as np
import os
import subprocess

def run(base_dir:str=os.path.dirname(os.path.abspath(__file__)), hdf_file:str=None, src_dataset:str=None, proj=None):
    cmd:str                 = None
    abs_hdf_file:str        = os.path.join(base_dir, hdf_file)
    name, exit              = os.path.splitext(hdf_file)

    geotiff_file:str        = "{0}_temp.tif".format(name)
    abs_geotiff_file:str    = os.path.join(base_dir, geotiff_file)
    
    dataset                 = gdal.Open(abs_hdf_file, gdal.GA_ReadOnly)
    metadata                = dataset.GetMetadata()
    geo_longs:list          = metadata['GRINGPOINTLONGITUDE.1'].split(', ')
    geo_lats:list           = metadata['GRINGPOINTLATITUDE.1'].split(', ')

    upper_long:str        = geo_longs[0]
    upper_lat:str         = geo_lats[0]
    lower_long:str        = geo_longs[3]
    lower_lat:str         = geo_lats[3]
    type:str              = "Int16"  
    if src_dataset is None: src_dataset:str = dataset.GetSubDatasets()[5][0] 

    cmd = "gdal_translate -a_srs EPSG:4326 -ot {type} -b 1 -a_ullr {ulx} {uly} {rlx} {rly} -of GTiff {src_dataset} {dst_dataset} -a_nodata {value}".format(type=type, value=np.nan, ulx=upper_long, uly=upper_lat, rlx=lower_long, rly=lower_lat, src_dataset=src_dataset, dst_dataset=abs_geotiff_file)

    result = subprocess.run(cmd.split(), capture_output=True)            

    dataset = None

if __name__=="__main__" :
    run(hdf_file="MAM35S0.A2008010.0000.002.2008010173218.hdf", src_dataset='HDF4_EOS:EOS_SWATH:"MAM35S0.A2008010.0000.002.2008010173218.hdf":mod35:Cloud_Mask')

This is the result of uploading the GTiff file converted using the above code to cesium ion.

However, I was surprised to see the results by uploading the results to cesium ion.

The image was projected inverted left and right. If anyone has experience or knowledge about this phenomenon, please help!

I also considered adding gcp information by duplicating the gdal_translate command.

But it didn't work, I suspected the projection was the problem, so I tried several different projections, but still no success.

I want to convert MODIS sensor's HDF file to GeoTiff so that it can be normally projected on GIS like cesium ion.

Solution

  • The swath data is ungridded, you cannot use gdal_translate for that. You can either use the GCP's in the metadata (eg with gdalwarp) or download the accompanying geolocation data and use that.

    It's not just that your data is flipped horizontally, look for example at the top near Crimea, that's a large shift as well. The use of -a_ullr overrides whatever the actual geolocation is, so it should be used with care.

    Using the GCP's with gdalwarp can for example be done with:

    gdalwarp -t_srs EPSG:4326 -tr 0.01 0.01 -te 25.8 28.4 33.7 46.9 HDF4_EOS:EOS_SWATH:"MAM35S0.A2008010.0000.002.2008010173218.hdf":mod35:Cloud_Mask cloudmask.tif

    In the end you'll have to decide what the output/target grid looks like, this is just a guess.

    If you're going to call it from Python, I'd recommend using the Python interface, so gdal.Translate or gdal.Warp etc. But that shouldn't influence the result.

    Here's a quicklook of the output overlayed on Google Maps using QGIS:

    enter image description here