[GIS] Convert MODIS mod021km hdf file to geotif using specific band

convertgdalmodispymodispython

I need to convert some bands of "mod021km" or "myd021km" files to geotif images (bands like: red,green, blue, nir ,…).

I used pymodis ==> modis_convert.py script but there is a problem!

When i select a subdataset, for example "EV_500_Aggr1km_RefSB" that it has 5 bands, as a result, I get a .TIF file with just one band!!! When I need a tif file for each band or one file with 5 bands.
I found out that pymodis uses gdal.ReprojectImage to convert image from HDF to Geotiff. i also found the code below in order to reproject hdf files:

from osgeo import gdal
import numpy as np

def hdf_subdataset_extraction(hdf_file, dst_dir, subdataset):
    """unpack a single subdataset from a HDF5 container and write to GeoTiff"""
    # open the dataset
    hdf_ds = gdal.Open(hdf_file, gdal.GA_ReadOnly)
    band_ds = gdal.Open(hdf_ds.GetSubDatasets()[subdataset][0], gdal.GA_ReadOnly)

    # read into numpy array
    band_array = band_ds.ReadAsArray().astype(np.int16)

    # convert no_data values
    band_array[band_array == -28672] = -32768

    # build output path
    band_path = os.path.join(dst_dir, os.path.basename(os.path.splitext(hdf_file)[0]) + 
    "-sd" + str(subdataset+1) + ".tif")

    # write raster
    out_ds = gdal.GetDriverByName('GTiff').Create(band_path,
                                              band_ds.RasterXSize,
                                              band_ds.RasterYSize,
                                              1,  #Number of bands
                                              gdal.GDT_Int16,
                                              ['COMPRESS=LZW', 'TILED=YES'])
    out_ds.SetGeoTransform(band_ds.GetGeoTransform())
    out_ds.SetProjection(band_ds.GetProjection())
    out_ds.GetRasterBand(1).WriteArray(band_array)
    out_ds.GetRasterBand(1).SetNoDataValue(-32768)

    out_ds = None  #close dataset to write to disc

One of the parameters of defining out_ds by gdal.GetDriverByName. Create is 1 (described as #number of bands). i changed it from 1 to 5 (number of subdataset bands) in modis_convert.py script and i get a .TIF file with 5 bands as a result. but i don't know that everything is correctly done or not! can anyone help me please.

Best Answer

im not familiar with MOD021 but you can try the code below. It will convert every band in the product. You need to loop across the hdf layers. If the subdataset is still a multiband then you need to loop in that particular dataset and use GetRasterBand. I have tried it and it works. This code will transform each band to TIFF. The names will be 'band1.tif', 'band2.tif' etc. Hope it helps.

from osgeo import gdal
import os


folder = 'C:/Users/Javier/Pecaries'
hdf_file = 'C:/Users/Javier/Pecaries/MOD021KM.A2020114.1445.061.2020115011536.hdf'

# open the dataset
hdf_ds = gdal.Open(hdf_file, gdal.GA_ReadOnly)
subdatasets = hdf_ds.GetSubDatasets()

for i in range(0, len(subdatasets)+1):
    subdataset_name = subdatasets[i][0]
    band_ds = gdal.Open(subdataset_name, gdal.GA_ReadOnly)
    band_path = os.path.join(folder, 'band{0}.tif'.format(i))
    if band_ds.RasterCount > 1:
        for j in range(1,band_ds.RasterCount + 1):
            band = band_ds.GetRasterBand(j)
            band_array = band.ReadAsArray()

    else:
        band_array = band_ds.ReadAsArray()
    out_ds = gdal.GetDriverByName('GTiff').Create(band_path,
                                                  band_ds.RasterXSize,
                                                  band_ds.RasterYSize,
                                                  1,
                                                  gdal.GDT_Int16,
                                                  ['COMPRESS=LZW', 'TILED=YES'])
    out_ds.SetGeoTransform(band_ds.GetGeoTransform())
    out_ds.SetProjection(band_ds.GetProjection())
    out_ds.GetRasterBand(1).WriteArray(band_array)
    out_ds.GetRasterBand(1).SetNoDataValue(-32768)



out_ds = None  #close dataset to write to disc
Related Question