Computación en la media de todos los rásteres en un directorio usando python

Question

No tengo ArcGIS (ArcPy), así que estoy esperando a usar algún otro paquete de python como numpy para calcular la media de una lista de los rásteres en un directorio (es decir, Raster1.tif + Raster2.tif + Raster3.tif / 3).

Me gusta la entrada de directorio en la consola para sus cálculos - el uso de este:

import os
import sys

directory = sys.argv[1]

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

Después de esto no estoy muy seguro de qué paquete o el código a utilizar. Los únicos ejemplos que he encontrado en internet son el uso de ArcPy.

Alguien me puede ayudar?

Yo estoy esperando a la exportación de un único promedio de trama como de 16 bits tif.

Answer 1

Usted podría utilizar rasterio leer sus datos como arrays de numpy (y escribir a partir de arrays de numpy) y numpy para realizar el cálculo del promedio.

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)

Answer 2

Para calcular la media de una lista de los rásteres en un directorio (es decir, Raster1.tif + Raster2.tif + Raster3.tif / 3), podría ser utilizado sólo una parte de Loïc Dutrieux rasterio código. Por ejemplo, en mi '/home/zeito/pyqgis_data' del directorio, el siguiente código:

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) 

produce en la Consola de Python de QGIS el siguiente resultado:

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

Si tienes problemas para instalar rasterio en su sistema, usted probablemente necesitará GDAL. El siguiente código produce el mismo resultado que en el código anterior:

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)