Creates a mosaic from some MODIS tiles for all SubDataset for a specific product

mosaic.py

# Creates a mosaic for the UK for a defined bounding box:
# - for all Days of Year (DoY) in julian day format for a specific year
# - using the tiles specified
# - for the selected product and each associated subdataset

import os
import gdal
import subprocess
from glob import glob

# This variables should come from the product yaml
# ================================================
tiles = ['h17v02', 'h17v03', 'h17v04', 'h18v03']
product = 'MOD13A2'
version = '006'
dataDir = '/data/peatlands_test_data/MOD13A2'
year = 2019

# Target extent
xmin = -784315
ymin = 5538675
xmax = 149109
ymax = 6786719
# ================================================

# ========================================================================
# From src/datacube/utils.py
def run_command(cmd: str):
    """
    Executes a command in the OS shell
    :param cmd: command to execute
    :return N/A:
    :raise Exception: If command executions fails
    """
    status = subprocess.call([cmd], shell=True)

    if status != 0:
        err_msg = f"{cmd} \n Failed"

def get_gdalbuildvrt_path():
    """
    Method to return the conda path for the gdalbuildvrt bin given the
    activated environment for the current python run.
    :return: Full path of gdalbuildvrt bin
    """

    conda_prefix = os.environ['CONDA_PREFIX']
    gdal_bin_path = os.path.join(conda_prefix, 'bin')
    return os.path.join(gdal_bin_path, 'gdalbuildvrt')
# ========================================================================

def get_number_of_subdatasets(fname):
    """
    Get number of available SubDatasets
    """
    d = gdal.Open(fname)
    sd = d.GetSubDatasets()

    return len(sd)

def get_subdataset(fname, index):
    """
    Get especific SubDataset by index
    """
    d = gdal.Open(fname)
    sd = d.GetSubDatasets()[index]

    return sd[0]

def get_filename(DoY, tile):
    """
    Get the full path file name corresponding to a DoY and tile
    """
    _fname = f'{product}.A{year}{DoY:03}.{tile}.{version}.*.hdf'
    _fname = glob(os.path.join(dataDir, _fname))

    if len(_fname) > 0:
        return _fname[0]
    else:
        return None

def get_sds_name(sds):
    """
    Get SubDataset name, replace spaces by underscores when needed,
    and remove any double quotes as well
    """
    return sds.split(':')[-1].replace(' ', '_').replace('"', '')

# Loop through all days of year (DoY) to get files and perform the mosaic
for DoY in range(1, 367):
    fnames = []
    for tile in tiles:
        _fname = get_filename(DoY, tile)
        if _fname is not None:
            fnames.append(_fname)

    # Check that we have one filename per tile
    if len(fnames) != len(tiles):
        # TODO Raise an error/warning...
        print(f'Missing files for DoY: {DoY:03}')
        continue

    # Get number of available SubDatasets
    n_sds = get_number_of_subdatasets(fnames[0])

    # Create mosaic for each SubDataset
    for sd in range(n_sds):
        # Get all SubDatasets
        sds = [f'"{get_subdataset(fname, sd)}"' for fname in fnames]
        # Get SDS name from first SDS
        sds_name = get_sds_name(sds[0])
        # List to string separated by spaces
        sds = ' '.join(sds)

        # Set outoput filename
        output_fname = f'{product}.A{year}{DoY:03}.{version}.{sds_name}.vrt'
        output_fname = os.path.join(dataDir, output_fname)

        command = (f'{get_gdalbuildvrt_path()} '
                   f'-te {xmin} {ymin} {xmax} {ymax} '
                   f'-overwrite {output_fname} {sds}')

        run_command(command)