Stacking a raster with GDAL in python - vrt, translate, merge giving errors

Question

I'm trying to create a composite of all my sentinel bands into one tif in python. I've tried doing a few things but nothing seems to be working. I am using an arcpy python environment as I need a few ArcGIS-specific tools later on in my script.

I have a dictionary of bands that are the band name and the absolute path to the tif.

First I tried using subprocess.call and gdal_merge.py

def createindices(index_dict, name):

    print(index_dict)

    coastal = index_dict["COASTAL_AEROSOL"]
    blue = index_dict["BLUE"]
    green = index_dict["GREEN"]
    red = index_dict["RED"]
    vegre5 = index_dict["VEG_RED_EDGE5"]
    vegre6 = index_dict["VEG_RED_EDGE6"]
    vegre7 = index_dict["VEG_RED_EDGE7"]
    nir = index_dict["NIR"]
    narrownir = index_dict["NARROW_NIR"]
    watervap = index_dict["WATER_VAPOUR"]
    swir11 = index_dict["SWIR11"]
    swir12 = index_dict["SWIR12"]

# Create the stacked image that has the highest resolution
    out = name + "_ALL_BANDS.tif"
    vrt = name + "_vrt.vrt"
    bandslist = [coastal, blue, green, red, vegre5, vegre6, vegre7, nir, narrownir, watervap, swir11, swir12]

    print("Stacking all bands")
    subprocess.call("gal_merge.py -o " + out + " -of GTiff -ps 10 10 -separate -a_nodata -10000 " + coastal + " " + blue + " " + green + " " +
                 red + " " + vegre5 + " " + vegre6 + " " + vegre7 + " " + nir + " " + narrownir + " " + watervap + " " + swir11
                + " " + swir12)

I get this error: FileNotFoundError: [WinError 2] The system cannot find the file specified. I am using the exact paths for my inputs, have checked it multiple times and they're right. I did a little reading online about subprocess giving these errors but I don't really understand the solutions.

I also tried replacing subprocees.call with os.system and get this error: 'gal_merge.py' is not recognized as an internal or external command, operable program or batch file.

So then I tried a different way to avoid using subprocess by using creating a vrt and then use gdal translate from this thread:

def createindices(index_dict, name):

    print(index_dict)

    coastal = index_dict["COASTAL_AEROSOL"]
    blue = index_dict["BLUE"]
    green = index_dict["GREEN"]
    red = index_dict["RED"]
    vegre5 = index_dict["VEG_RED_EDGE5"]
    vegre6 = index_dict["VEG_RED_EDGE6"]
    vegre7 = index_dict["VEG_RED_EDGE7"]
    nir = index_dict["NIR"]
    narrownir = index_dict["NARROW_NIR"]
    watervap = index_dict["WATER_VAPOUR"]
    swir11 = index_dict["SWIR11"]
    swir12 = index_dict["SWIR12"]

# Create the stacked image that has the highest resolution
    out = name + "_ALL_BANDS.tif"
    vrt = name + "_vrt.vrt"
    bandslist = [coastal, blue, green, red, vegre5, vegre6, vegre7, nir, narrownir, watervap, swir11, swir12]

    print("Stacking all bands")

    vrt_options = gdal.BuildVRTOptions(resolution='10', separate=True)
    gdal.BuildVRT(vrt, bandslist, options=vrt_options)

    translate_options = gdal.TranslateOptions(format="GTiff", outputType='gdal.GDT_Int16', noData=-'10000')
    gdal.Translate(out, vrt, options=translate_options)

I get this error: TypeError: in method 'BuildVRTInternalNames', argument 3 of type 'GDALBuildVRTOptions *'

I also can't use rasterio as I can't duplicate my arcpy environment for some reason, which you need to do to add new libraries to it.

I'm lost.

Answer 1

Could also be done with:

gdalbuildvrt -input_file_list ListOfFilesToInclude.list -separate NameOfVrt.vrt

Answer 2

Managed to solve this as I was typing out my question, so I decided to share my work around in case anyone else has this issue:

I tried this from looking at this site:

def createindices(index_dict, name):

    print(index_dict)

    coastal = index_dict["COASTAL_AEROSOL"]
    blue = index_dict["BLUE"]
    green = index_dict["GREEN"]
    red = index_dict["RED"]
    vegre5 = index_dict["VEG_RED_EDGE5"]
    vegre6 = index_dict["VEG_RED_EDGE6"]
    vegre7 = index_dict["VEG_RED_EDGE7"]
    nir = index_dict["NIR"]
    narrownir = index_dict["NARROW_NIR"]
    watervap = index_dict["WATER_VAPOUR"]
    swir11 = index_dict["SWIR11"]
    swir12 = index_dict["SWIR12"]

# Create the stacked image that has the highest resolution
    out = name + "_ALL_BANDS.tif"
    vrt = name + "_vrt.vrt"
    bandslist = [coastal, blue, green, red, vegre5, vegre6, vegre7, nir, narrownir, watervap, swir11, swir12]

    print("Stacking all bands")
# make vrt
    destName = vrt
    srcDSOrSrcDSTab = bandslist

    kwargs = {
        'xRes': '10',
        'yRes': '10',
        'separate': True
    }

    ds = gdal.BuildVRT(destName, srcDSOrSrcDSTab, **kwargs)

# close and save ds
    ds = None

# save vrt to tif with gdal translate
    kwargs = {
        'format': 'GTiff',
        'outputType': gdal.GDT_Int16,
        'noData': '-10000'
    }

    fn = vrt
    dst_fn = out

    ds = gdal.Translate(dst_fn, fn, **kwargs)

# close and save ds
    ds = None