Note that the Command Line method only adds the main images to the multiband GTiff without the H5s additional metadata. The Python method is more complicated, but it will place all subdatasets into the GTiff, as well as the additional related metadata.
Using the Command Line
file="path/to/input.h5"
subdataset1="${file}://science/LSAR/GCOV/grids/frequencyA/HHHH"
subdataset2="${file}://science/LSAR/GCOV/grids/frequencyA/numberOfLooks"
subdataset3="${file}://science/LSAR/GCOV/grids/frequencyA/rtcAreaNormalizationFactorGamma0ToSigma0"
dimensions_args="-ul_lr upper_left_corner_x upper_left_corner_y lower_right_corner_x lower_right_corner_y -ps pixel_size_x pixel_size_y"
gdal_merge.py $dimensions_args -seperate -of GTiff -o $(echo $file | sed 's/\.h5/\_multiband.tif/g')
Using Python
Here is the code to extract raster subsets in h5 and convert them into geotiff files.
import json
import subprocess
import re
from osgeo import gdal
file="NISAR_L2_PR_GCOV_001_005_A_219_4020_SHNA_A_20081012T060910_20081012T060926_P01101_F_N_J_001.h5"
actual_json = json.loads(subprocess.check_output(f"gdalinfo -json {file}", shell=True))
datasets = actual_json['metadata']['SUBDATASETS']
names = [key for key in datasets.keys() if re.match('SUBDATASET_[0-9]*_NAME', key)]
for name in names:
subdata = f"{datasets[name].replace('HDF5','NETCDF')}"
if gdal.Open(subdata).GetProjectionRef() != '':
outfile =f"{subdata.replace('/','_')}.tif"
subprocess.run(f'gdal_translate -of GTiff {subdata} {outfile}', shell=True)
We also provide the code to combine all raster subdatasets with the same datatype, extent, size, and projection in the h5 into multiple-band geotiff files. We notice that the combined multiple band geotiff files do not include the metadata for each band in the geotiff files.
import json
import subprocess
import re
import numpy as np
from osgeo import gdal
def combine_geotiff(filelist, outfile):
'''
Merge the geotiff files with the same data type, extent, and projection into a multiple band geotif file
:param filelist: list of geotiff files
:param outfile: filename of the output file
:return: None
'''
bandmetalist = [gdal.Info(file, format='json')['bands'] for file in filelist]
bandnumlist = [len(item) for item in bandmetalist]
bandnum = int(np.array(bandnumlist).sum())
filenum = len(filelist)
for i in range(filenum):
ds = gdal.Open(filelist[i])
if i == 0:
# create output raster file
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(outfile, ds.RasterXSize, ds.RasterYSize, bandnum, ds.GetRasterBand(1).DataType )
outdata.SetGeoTransform(ds.GetGeoTransform()) ##sets same geotransform as input
outdata.SetProjection(ds.GetProjection()) ##sets same projection as input
bandnum_of_file = ds.RasterCount
fmeta = ds.GetMetadata()
for j in range(1, bandnum_of_file + 1):
band = ds.GetRasterBand(j)
outdata.GetRasterBand(i*bandnum_of_file + j).WriteArray(band.ReadAsArray())
fmeta.update(band.GetMetadata())
outdata.GetRasterBand(i*bandnum_of_file + j).SetMetadata(fmeta)
ds = None
outdata = None
file="NISAR_L2_PR_GCOV_001_005_A_219_4020_SHNA_A_20081012T060910_20081012T060926_P01101_F_N_J_001.h5"
actual_json = json.loads(subprocess.check_output(f"gdalinfo -json {file}", shell=True))
datasets = actual_json['metadata']['SUBDATASETS']
names = [key for key in datasets.keys() if re.match('SUBDATASET_[0-9]*_NAME', key)]
datalist=[]
for name in names:
subdata = f"{datasets[name].replace('HDF5','NETCDF')}"
ds = gdal.Open(subdata)
if ds.GetProjectionRef() != '':
size = f'{ds.RasterXSize}x{ds.RasterYSize}'
outfile = f"{subdata.replace('/','_')}.tif"
outfile = outfile.split(":")[2].replace('__','')
mo=outfile.replace(".tif","")
subprocess.run(f"gdal_translate -of GTiff -ot Float32 -mo Band={mo} {subdata} {outfile}", shell=True)
datalist.append({"filename":outfile, "size":size})
sizelist = [item['size'] for item in datalist]
uniqvals= np.unique(sizelist)
for item in uniqvals:
itemlist = [ dataitem['filename'] for dataitem in datalist if dataitem['size'] == item ]
mb_outfile = subdata.split(":")[1].replace('.h5',f'_{item}.tif')
mb_outfile = mb_outfile.replace('"','')
combine_geotiff(itemlist, mb_outfile)
