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.