[GIS] Computing mean of all rasters in a directory using python

pythonraster

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.

Best Answer

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)
Related Question