<g id="1">¿Cuá</g><g id="2">l es la mejor forma de sacarle el m<g id="3">áximo partido a las tendencias?</g></g>

Question

Estoy buscando el enfoque más sencillo tomar la media a través de múltiples imágenes rasterizadas, que contienen una gran cantidad de valores de nan. El valor promedio debe ignorar cualquier nan. El gdal_calc no funciona, porque salidas nodata cuando cualquier valor nodata es encoutered.

¿Hay un enfoque simple con GDAL o python (con cualquier paquete)? Prefiero no usar ArcGIS.

Answer 1

¿Usted está mirando para tomar el promedio para cada celda de la cuadrícula en la pila, o el promedio general? Si es la primera, se podría utilizar el AverageOverlay herramienta en la WhiteboxTools de la biblioteca. Esto puede ser escrito en Python como sigue:

from whitebox_tools import WhiteboxTools

wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"

wbt.average_overlay(inputs='file1.tif;file2.tif;file3.tif', output='average.tif')

Todas las herramientas en esta biblioteca se manejar valores NoData adecuadamente suponiendo que el valor NoData está marcado correctamente en los archivos. Para obtener más información sobre el AverageOverlay herramienta, consulte la WhiteboxTools Manual de Usuario. Para la herramienta de código fuente, ver aquí, y en particular, observe cómo NoData valores que se manejan en esta línea en el código. Es decir, la herramienta calcula un promedio para cada celda de la cuadrícula sólo se basa en los valores válidos en la pila de imágenes. El WhiteboxTools binario ejecutable puede ser descargado desde el Geomorphometry y Hydrogeomatics Grupo de Investigación de la página. Hay un WhiteboxTools plugin de QGIS, aunque aún es algo experimental y se actualizará en breve.

Answer 2

Ya se indicó que la lectura de los rásteres en la memoria estaría bien, he aquí un ejemplo simple de cómo puede obtener su promedio con gdal y numpy. Supongo que todos los que aquí su tiff son del mismo tamaño (filas, columnas), son solo bandas, comparten el mismo CRS, y usted no necesita escribir la calculada imagen de disco.

import gdal
import numpy as np

def write_raster(raster_array, gt, data_obj, outputpath, dtype=gdal.GDT_UInt16, options=0, color_table=0, nbands=1, nodata=False):

    height, width = raster_array.shape

    # Prepare destination file
    driver = gdal.GetDriverByName("GTiff")
    if options != 0:
        dest = driver.Create(outputpath, width, height, nbands, dtype, options)
    else:
        dest = driver.Create(outputpath, width, height, nbands, dtype)

    # Write output raster
    if color_table != 0:
        dest.GetRasterBand(1).SetColorTable(color_table)

    dest.GetRasterBand(1).WriteArray(raster_array)

    if nodata is not False:
        dest.GetRasterBand(1).SetNoDataValue(nodata)

    # Set transform and projection
    dest.SetGeoTransform(gt)
    wkt = data_obj.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    dest.SetProjection(srs.ExportToWkt())

    # Close output raster dataset 
    dest = None

tifflist = ['something1.tif', 'something2.tif'] # the paths to each of the tiffs you want to use in averaging
for i, tiff in enumerate(tifflist):
    gd_obj = gdal.Open(tiff)
    array = gd_obj.ReadAsArray()
    array = np.expand_dims(array,2)
    if i == 0:
        allarrays = array
    else:
        allarrays = np.concatenate((allarrays, array), axis=2)
mean_of_tiffs = np.nanmean(allarrays, axis=2)

outputpath = 'wherever you are saving this guy.tif'
write_raster(mean_of_tiffs, gd_obj.GetGeoTransform(), gd_obj, outputpath)