19

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?

Aaron
  • 51,658
  • 28
  • 154
  • 317
Jens Hiestermann
  • 365
  • 1
  • 3
  • 5

4 Answers4

34

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

Loïc Dutrieux
  • 2,302
  • 17
  • 17
  • 1
    Yes, this is essentially what Rasterio's rio-stack program does: https://github.com/mapbox/rasterio/blob/master/rasterio/rio/stack.py#L99-L110. – sgillies Jan 11 '17 at 07:25
  • 1
    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 '19 at 10:21
  • alas I get this error "RasterioIOError: '/' not recognized as a supported file format." – ilFonta Jul 31 '19 at 19:37
  • @ilFonta, make a new question with a reproducible, minimal code example then. – bugmenot123 Jan 08 '20 at 15:41
  • @sgillies I tried a few passes at using rio stack but couldn't figure out how to get it to work in this example. Assuming we have 3 tifs R,G, and B, how would it work? I'm getting squat – mmann1123 Jan 25 '22 at 01:45
16

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)
user2856
  • 65,736
  • 6
  • 115
  • 196
  • if the every file of tifs list is mutli-bands file, and diferent band in every files need to be stacked to outfile, meanwhile sparate is True. This code is not work – user7177639 Jun 14 '20 at 08:29
  • 1
    That scenario is not part of the question. If you have a new question, then ask it separately. – user2856 Jun 14 '20 at 10:48
  • i am sorry to rashly comment. my purpose is to extend this question for solving these kind of layer stacking by gdal. I am sorry again – user7177639 Jun 14 '20 at 11:17
0

With gdal and numpy the following is possible assuming the images are the same size (rows and columns) and have the same CRS (e.g. 10m bands from a Sentinel 2 scene).

import gdal
import numpy as np

file_list = ['band1.tif', 'band2.tif', 'band3.tif']


array_list = []


Read arrays


for file in file_list:
    src = gdal.Open(file)
    geotransform = src.GetGeoTransform() # Could be done more elegantly outside the for loop
    projection = src.GetProjectionRef()
    array_list.append(src.ReadAsArray())
    src = None


Stack arrays


stacked_array = np.stack(array_list, axis=0)
array_list = None


Write to disk


driver = gdal.GetDriverByName('GTiff')
n, rows, cols = stacked_array.shape
dataset = driver.Create('output_file_name.tif', cols, rows, n,
                                gdal.GDT_Uint16)
dataset.SetGeoTransform(geotransform)
dataset.SetProjection(projection)


for b in range(1,n+1):
    band = dataset.GetRasterBand(b) # GetRasterBand is not zero indexed
    band.WriteArray(stacked_array[b-1]) # Numpy is zero indexed


dataset = None
stacked_array = None
PyMapr
  • 1,630
  • 1
  • 16
  • 31
0

(thanks to https://gis.stackexchange.com/a/299218/155670, https://gis.stackexchange.com/a/223915/155670)

import osgeo.gdal
from osgeo.gdal import Dataset

vrt_internal_name: str = '' # if non-empty and not the same as target_tif, vrt will be created
vrt: Dataset = osgeo.gdal.BuildVRT(vrt_internal_name, ['r_0.tif', 'r_1.tif', 'r_2.tif'])


target_tif: str = 'r.tif'
osgeo.gdal.Translate(target_tif, vrt)

For me to make BuildVRT work all paths (inputs and outputs) were needed to be resolved

Георгий Савельев
  • 311
  • 3
  • 5