5

I don't have ArcGIS (ArcPy) so I am hoping to use some other python package such as numpy to compute the mean of a list of rasters in a directory (i.e. Raster1.tif + Raster2.tif + Raster3.tif / 3).

I like to input the directory in the console for computations- using this:

import os
import sys

directory = sys.argv[1]

rasters []
for i in os.listdir(directory):
    if i.endswith(".tif"):

After this I am not really sure what package or code to use. The only examples I have found online are using ArcPy.

Can anyone help me out?

I am hoping to export a single averaged raster as a 16-bit tif.

  • Are the rasters single band images or multi-band? – Kersten Jun 18 '17 at 12:05
7

You could use rasterio to read your data as numpy arrays (and write from numpy arrays) and numpy to perform the averaging.

import rasterio
import numpy as np
from glob import glob
import os

data_dir = '/path/to/data' # Or sys.argv[1]
file_list = glob(os.path.join(data_dir, '*.tif'))

def read_file(file):
    with rasterio.open(file) as src:
        return(src.read(1))

# Read all data as a list of numpy arrays 
array_list = [read_file(x) for x in file_list]
# Perform averaging
array_out = np.mean(array_list, axis=0)

# Get metadata from one of the input files
with rasterio.open(file_list[0]) as src:
    meta = src.meta

meta.update(dtype=rasterio.float32)

# Write output file
with rasterio.open('/path/to/output/file.tif', 'w', **meta) as dst:
    dst.write(array_out.astype(rasterio.float32), 1)
  • Thanks Loic- this code is now very close to working. Not sure if you know what the source of my traceback error is... I used your code (small typo on line 21 fixed ('scr.meta'). data_dir was assigned sys.argv[1] and put into the output directory. here is my error – Alexander Jun 18 '17 at 16:29
  • C:\Users\Alex\Anaconda2\envs\alenv\lib\site-packages\rasterio_init_.py:282: NotGeoreferencedWarning: Dataset has no geotransform set. Default transform will be applied (Affine.identity()) s.start() Traceback (most recent call last): File "data_processing/flatfield_compute.py", line 26, in <module> with rasterio.open(data_dir, 'w', **meta) as dst: File "C:\Users\Alex\Anaconda2\envs\alexenv\lib\site-packages\rasterio_init_.py", line 282, in open s.start() File "rasterio_io.pyx", line 986, in rasterio._io.DatasetWriterBase.start (rasterio/_io.c:17228) – Alexander Jun 18 '17 at 16:34
  • Thanks for the typo. Is that really an error you're getting or 'just' a warning? It looks like your input files aren't geo-referenced. – Loïc Dutrieux Jun 19 '17 at 16:04
  • Hey Loic- yes. My files aren't georeferenced because I am computing flatfields with this script. So my data do not need spatial info. I had tried hashtagging out the metadata part of the script but it still doesn't seem to run. I'll keep playing around with it unless you have some suggestions. – Alexander Jun 21 '17 at 11:56
  • I just tried with simple pictures taken with my phone (not georeferenced), and the process completes without any problem. Are you sure what you are getting is an error and not a warning? Do all your image have the same exact dimensions? – Loïc Dutrieux Jun 21 '17 at 16:22
3

To compute the mean of a list of rasters in a directory (i.e. Raster1.tif + Raster2.tif + Raster3.tif / 3), it could be used only a part of Loïc Dutrieux rasterio code. For example, in my '/home/zeito/pyqgis_data' directory, next code:

import rasterio
import numpy as np
from glob import glob
import os

data_dir = '/home/zeito/pyqgis_data' # Or sys.argv[1]
file_list = glob(os.path.join(data_dir, '*.tif'))
raster = [ os.path.split(item)[1] for item in file_list ]

def read_file(file):
    with rasterio.open(file) as src:
        return(src.read(1))

# Read all data as a list of numpy arrays 
array_list = [read_file(x) for x in file_list]
# Perform averaging
for i, array in enumerate(array_list):
    print raster[i], ", mean: ", np.mean(array) 

produces at the Python Console of QGIS the following result:

utah_lake.tif , mean:  116.720311851
b4.ND.tif , mean:  68.2443939062
LT50380322011235PAC01_B6.tif , mean:  162.756999634
b3.ND.tif , mean:  42.4123148545
tiznados_canoa_part_reproj.tif , mean:  610.861699623
LT50380322011235PAC01_B3.tif , mean:  43.0117399461
tiznados_canoa_part.tif , mean:  633.044776
natural_earth.tif , mean:  152.838815689
LT50380322011235PAC01_B4.tif , mean:  68.3190660858
tiznados_canoa.tif , mean:  256.409105513
utah_demUTM2.tif , mean:  1824.71800614

If you have issues to install rasterio in your system, you probably will need GDAL. Next code produces the same result as above code:

from osgeo import gdal
import numpy as np
from glob import glob
import os

data_dir = '/home/zeito/pyqgis_data' # Or sys.argv[1]
file_list = glob(os.path.join(data_dir, '*.tif'))
raster = [ os.path.split(item)[1] for item in file_list ]

for i, file in enumerate(file_list):
    dataset = gdal.Open(file)
    band = dataset.GetRasterBand(1)
    data = band.ReadAsArray()

    print raster[i], ", mean: ", np.mean(data)
  • thanks xunilk! This code also works! I prefer to output my file as a tiff rather than printing the mean values. However, this code will also be very useful to me for verifying the means before averaging. Thanks so much \ – Alexander Jun 18 '17 at 16:35

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.