¿Agregar raster con GDAL?

Question

Tengo un GeoTIFF con una resolución de 30 metros. Necesito reducir la muestra a una resolución de 250m (utilizando alguna función de agregación como media, modo, suma). Mi pila es python (rasterio y GDAL bindings)

Puedo encontrar una manera de hacer esto directamente usando numpy escribiendo algunos bucles dobles desagradables e iterando sobre las matrices directamente, pero buscando una manera más limpia.

Tenga en cuenta que el principal problema es que 250 no es un múltiplo de 30 y habrá algún desbordamiento. Por ejemplo si la función de agregación es la suma el algoritmo sería utilizar:

1. Bloque de 8 X 8 que tiene una resolución base de 240 m X 240 m
2. Elige una franja parcial de píxeles a 10m X 10m para completar la suma a 250m
3. El valor de 2. se ponderaría por área y se añadiría a la suma de 1.

Por favor, indíqueme si no he sido claro. Una vez más, puedo encontrar una forma de fuerza bruta para hacerlo utilizando manipulaciones de la matriz, pero me preguntaba si hay una manera elegante. Estoy abierto a usar R también y leer sobre el projectRaster y los métodos agregados, pero parece que también funcionan con cortes enteros.

En general, puede ser realmente engorroso hacer esto para una resolución arbitraria, por ejemplo, 31,567 m, más o menos.

Answer 1

Yo usaría gdalwarp para remuestrear (reducir la muestra en este caso) los datos:

gdalwarp -tr 250 250 -r average input_r30.tif output_r250.tif

En su lugar, utiliza los enlaces de GDAL Python:

from osgeo import gdal

inDs = gdal.Open('input_r30.tif')
outDs = gdal.Warp('output_r250.tif', inDs,
                  format = 'GTiff',
                  xRes = 250, yRes = 250,
                  resampleAlg = gdal.GRA_Average)

inDs = None
outDs = None

En ambos ejemplos, he utilizado el método de remuestreo medio, pero puedes elegirlo entre los admitidos:

GRA_NearestNeighbour = 0
GRA_Bilinear = 1
GRA_Cubic = 2
GRA_CubicSpline = 3
GRA_Lanczos = 4
GRA_Average = 5
GRA_Mode = 6
GRA_Gauss = 7 

Como 250 no es un múltiplo de 30, la salida tendrá una extensión espacial un poco mayor que la entrada.

---

ACTUALIZACIÓN: Gracias a @MikeT, que ha abierto un billete en el rastreador de errores de GDAL, "las constantes GRA adicionales estarán disponibles por conveniencia a través de los enlaces de SWIG Python en el próximo Lanzamiento de GDAL 2.2.3 pronto estará disponible". Estas son las constantes relativamente nuevas:

GRA_Max = 8
GRA_Min = 9
GRA_Med = 10
GRA_Q1 = 11
GRA_Q3 = 12

---

ACTUALIZACIÓN: GDAL 2.2.3 ha sido publicado hoy (2017/11/23):
https://lists.osgeo.org/pipermail/gdal-dev/2017-November/047782.html