0

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.

2 Answers 2

1

Could also be done with:

gdalbuildvrt -input_file_list ListOfFilesToInclude.list -separate NameOfVrt.vrt
0

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
1
  • 2
    could you clarify as to what the solution actually is?
    – Ian Turton
    Mar 6, 2020 at 14:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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