2

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?

2
  • Is the maximum number of rasters to stack 4 or it's variable?
    – xunilk
    Commented Mar 14, 2017 at 18:54
  • I see a problem in the order of bands. I'm afraid that the ordering of natsort was not good (seems that bands to stack are not adjacent). So, I need a complete series of your bands to stack for cracking the algorithm. Probably you need to use re (regular expressions) python module instead natsort.
    – xunilk
    Commented Mar 14, 2017 at 23:08

2 Answers 2

2

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.

enter image description here

10
  • Dear xunilk, thank you for your help, but unfortunately it doesn't work. Nothing happens. I don't know why
    – ilFonta
    Commented Mar 13, 2017 at 13:02
  • Obviously, its works with my criteria. You need to establish your own stacking criteria. To find out yours, please, comment all lines after for i in range(1, n, 2) (including this one), run the script and put in your question an extract of obtained ord_list.
    – xunilk
    Commented Mar 13, 2017 at 16:14
  • How many bands do you want to stack? In my case two bands was only an example.
    – xunilk
    Commented Mar 13, 2017 at 18:49
  • Dear Xunilk, your scripts works, but doesn't properly solve my problem. I don't want n stacks with two bands, but only one stack with n bands. I modified your script and i will report only the last part. n = len(ord_list) for i in range(0, n, 1): band = ord_list[i] a = sourceDir + band a.append(b) cmd = "gdal_merge.py -separate -o " + sourceDir + namestack + " " + ddd + " os.system(cmd) I receive an error a.append(b) AttributeError: 'str' object has no attribute 'append'
    – ilFonta
    Commented Mar 13, 2017 at 18:56
  • I understood but, I need an extract of your ord_list to see the structure of your file rasters.
    – xunilk
    Commented Mar 13, 2017 at 19:02
2

I've found the solution. Just this; @Xunilk

import os
import natsort

sourceDir = "/home/giacomo/Desktop/eee/"
unord_list = set(os.listdir(sourceDir))

ord_list = []
dd = ""
namestack = "prova.tif"

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(0, n, 1):
    band = ord_list[i]        
    aa = sourceDir + band
    dd = dd + aa + " " 

cmd = "gdal_merge.py -separate -o" + " "+ sourceDir + namestack + " " 
+ dd
os.system(cmd)

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.