[GIS] Raster stacking with naturally ordered layers gives ValueError: need at least one array to concatenate

Question

I need to create a multiple bands .tiff file, also known as rasters stack.
I ordered the files to be stacked (from the folder) using natsort and it works.

How to use the method 'os.walk' with the ordered list (ord_list), because
I get this error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/tmp/tmpAP_F30.py", line 20, in <module>
    stackRast = np.dstack(rasters)
  File "/usr/lib/python2.7/dist-packages/numpy/lib/shape_base.py", line 368, in dstack
    return _nx.concatenate([atleast_3d(_m) for _m in tup], 2)
ValueError: need at least one array to concatenate

I wrote my script doing copy/paste of some stuff in the net, I'm not skilled in Python yet.

import os
import natsort
import sys
import numpy as np
import gdal
import gdalnumeric as gd
from scipy.signal import savgol_filter

rasters = list()

# natural sorting
sourceDir = "/home/giacomo/Desktop/eee/"
unord_list = set(os.listdir(sourceDir))
ord_list = ('\n'.join(natsort.natsorted(unord_list)))

# layer stacking
for folder, subs, files in os.walk(ord_list):
    for filename in files:
        aSrc = gd.LoadFile(os.path.join(folder,filename))
        rasters.append(aSrc)

stackRast = np.dstack(rasters)

How can I avoid getting this error?

Answer 1

The issue is because sourceDir is a directory and unord_list is precisely a list (no a directory). Command os.walk(ord_list) cannot work in this context. You need to use ord_list only as criterion for matching the list of rasters to be stacked. By the way, for stacking your rasters you can use gdal_merge.py with -separate parameter. It's straightforward.

To test my approach, I modified your code to force my rasters b3_ref.tif and b4_ref.tif (in my sourceDir) to stack among 174 *tif files. You should modified it with your own criterion.

import os
import natsort

sourceDir = "/home/zeito/pyqgis_data/"
unord_list = set(os.listdir(sourceDir))
ord_list = []

for folder, subs, files in os.walk(sourceDir):
    for filename in files:
        base_name, extension = os.path.splitext(filename)
        if(extension == ".tif"):
            ord_list.append(filename)

ord_list = natsort.natsorted(ord_list)

n = len(ord_list)

for i in range(1, n, 2): #warning: iterator begins from 1 to match
                         #2 layers in my example
    try:
        b1 = ord_list[i] 
        b2 = ord_list[i+1]
        output = "stack_" + b1[:-4] + "_" + b2[:-4] + ".tif"
        if b1 == 'b3_ref.tif' and b2 == 'b4_ref.tif':
            cmd = 'gdal_merge.py -separate -o ' + sourceDir + output + ' ' + sourceDir + b1 + ' ' + sourceDir + b2
            os.system(cmd)
            print "done"
    except IndexError:
        pass

After running the code at the Python Console of QGIS, the stacked layer was named as stack_b3_ref_b4_ref.tif and it looks like as in next image. I used Value Tool plugin to corroborate that it's a stacked raster and its values are coincident with original rasters.