5

I'm trying to create a mosaic raster from two rasters. With rasterio works well but gdal_merge.py subprocess calling doesn`t work for me.

My Error is:

OSError: [WinError 193] %1 is not a valid Win32 application

import numpy as np
import matplotlib.pyplot as plt
import subprocess, glob
from osgeo import gdal

files_to_mosaic = glob.glob('C:/Users/DanielKovacs/Documents/GranCanaria/mosaic_L2A_T28RDS_A025431_20200505T115222_2020- 05-05_con/*_B04.tif')
files_to_mosaic

files_string = " ".join(files_to_mosaic)
print(files_string)

command = "gdal_merge.py -o mosaic2.tif -of gtiff " + files_string
output = subprocess.run(command)
output
6
  • subprocess cannot find gdal_merge.py, you have to specify full path to gdal_merge.py or use os.system which inherits the path from the os. You can erase imports for numpy and matplotlib and gdal.
    – Zoltan
    Commented May 10, 2020 at 13:58
  • The same happens command = "C:/ProgramData/Anaconda3/Scripts/gdal_merge.py -o TEAK_Aspect_Mosaic.tif -of gtiff " + files_string output = subprocess.run(command) output Commented May 10, 2020 at 13:59
  • command = "C:/ProgramData/Anaconda3/Scripts/gdal_merge.py -o C:/Users/DanielKovacs/Documents/GranCanaria/mosaic_L2A_T28RDS_A025431_20200505T115222_2020-05-05_con/output.tif -of gtiff " + files_string output = os.system(command) output With this statement gives me no error, however the output file is not created in my output path. Commented May 10, 2020 at 14:05
  • 3
    Why not just python gdal_merge.py -o out.tif in1.tif in2.tif as in the example gdal.org/programs/gdal_merge.html? Did you know that for most use cases (IMHO) gdalbuildvrt gdal.org/programs/gdalbuildvrt.html#gdalbuildvrt is a better option than gdal_merge?
    – user30184
    Commented May 10, 2020 at 18:47
  • 1
    Your're right! gdalbuildvrt is just amazing and fast. I had executed from OSGeo4W Shell Commented May 10, 2020 at 19:38

2 Answers 2

3

You can probably do this with gdal.Warp in Python directly:

from osgeo import gdal

files_to_mosaic = ["a.tif", "b.tif"] # However many you want.
g = gdal.Warp("output.tif", files_to_mosaic, format="GTiff",
              options=["COMPRESS=LZW", "TILED=YES"]) # if you want
g = None # Close file and flush to disk

There's loads of other options you can add on to this.

1
  • Very usefull, thank you! Commented May 18, 2020 at 13:31
0

I've used @Jose answer but I had issues with nodata (apparently, an option for telling gdal.Warp what values to ignore is necessary), I had the black of superior images covering the one below.

Using rasterio though worked nice (from this post: https://gis.stackexchange.com/a/247430/88415):

import rasterio
from rasterio.merge import merge

with rasterio.open(files_tif[0]) as src:
    meta = src.meta.copy()

# The merge function returns a single array and the affine transform info
arr, out_trans = merge(files_tif)

meta.update({
    "driver": "GTiff",
    "height": arr.shape[1],
    "width": arr.shape[2],
    "transform": out_trans
})

# Write the mosaic raster to disk
with rasterio.open(output, "w", **meta) as dest:
    dest.write(arr)

Where output = 'mosaic2.tif' and files_tif is the list of input tif files.

1
  • 1
    gdal.Warp has options srcNodata and dstNodata to e.g. translate all input missing values into (e.g.) a single output "no data" value.
    – Jose
    Commented May 21, 2021 at 13:30

Your Answer

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

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