3

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.

3
  • Hey, welcome to GIS SE! Do you have a file I can use to test? May 6, 2020 at 15:00
  • thanks @Marcelo for your response, for example "MOD021KM.A2019183.0705.061.2019183190722.hdf" is one of the files i'm trying to convert it's bands to geotiff. May 12, 2020 at 6:23
  • DO you have a link where I can download the file? May 13, 2020 at 13:35

1 Answer 1

3

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
3
  • thanks a lot @jfoguet. this works but i can't understand the difference between gdal.ReprojectImage method and the method that it is been used in the code above. which method should i use? May 11, 2020 at 11:06
  • hi @ Mehrdad HMT. Glad it worked. My response was a variation to the piece of code you found and tried to solve the problem of 'unpacking' every band in the HDF as a single band. The code also gets the HDF original geographic info (if there is one) and sets it to your new TIFF. (SetGeotransform and SetProjection, methods) gdal.ReprojectImage is usefull to reproject your data from one crs to another (in MODIS products typically from a sinusoidal projection to whatever projection you want to use)
    – jfoguet
    May 11, 2020 at 18:05
  • hi @jfoguet. with a little bit changes in your code, i could save every band as TIF file but there is another problem. since mod021km or myd021km files are not georef, so they don't have any geotransform or projection values. GetGeoTransform function returns incorrect values and GetProjection function returns empty string "". so i guess that the only way to georefrence modis bands is using gdal.ReprojectImage. but as far as i saied gdal.ReprojectImage does not work with single bands. it takes the layer with all of it's bands as a "dataset" format and reprojects it in a file with one band! May 13, 2020 at 11:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.