8

Does anyone have a eloquent way of stacking multiple .tif files into a multiple band stack using Rasterio and/or GDAL?

I am looking for a way to avoid using a subprocess command like gdal_merge.py and rather have it as part of my python script.

I know that both Rasterio and GDAL read the .tif files as arrays, but how do I stack those arrays and write out the result as separate bands?

14

Using rasterio you could do

import rasterio

file_list = ['file1.tif', 'file2.tif', 'file3.tif']

# Read metadata of first file
with rasterio.open(file_list[0]) as src0:
    meta = src0.meta

# Update meta to reflect the number of layers
meta.update(count = len(file_list))

# Read each layer and write it to stack
with rasterio.open('stack.tif', 'w', **meta) as dst:
    for id, layer in enumerate(file_list, start=1):
        with rasterio.open(layer) as src1:
            dst.write_band(id, src1.read(1))

It of course assumes that your input layers already share the same extent, resolution and data type

  • 1
    Yes, this is essentially what Rasterio's rio-stack program does: github.com/mapbox/rasterio/blob/master/rasterio/rio/…. – sgillies Jan 11 '17 at 7:25
  • Is it efficient to keep the stack in memory (to perform several functions on the various bands) instead of writing the stacked file out? Or should it be written to a file and then manipulated? – Shawn Jan 25 at 10:21
  • alas I get this error "RasterioIOError: '/' not recognized as a supported file format." – ilFonta Jul 31 at 19:37
8

If using GDAL 2.1+ it's as simple as gdal.BuildVRT then gdal.Translate:

from osgeo import gdal
outvrt = '/vsimem/stacked.vrt' #/vsimem is special in-memory virtual "directory"
outtif = '/tmp/stacked.tif'
tifs = ['a.tif', 'b.tif', 'c.tif', 'd.tif'] 
#or for all tifs in a dir
#import glob
#tifs = glob.glob('dir/*.tif')

outds = gdal.BuildVRT(outvrt, tifs, separate=True)
outds = gdal.Translate(outtif, outds)
  • It raises an error: ValueError: Received a NULL pointer and points to the Translate method – Rodrigo E. Principe May 30 '18 at 12:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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