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], 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)) + "-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.