Procesamiento de geotiff de gran tamaño con python

Question

Estoy trabajando con un geotiff grande (30GB) de 3 bandas. Estoy usando el siguiente código para intentar importar las bandas en python

import gdal
import numpy as np

raster_path = "C:/Path/To/Image.tif"
raster_dataset = gdal.OpenEx(raster_path, gdal.GA_ReadOnly)
geo_transform = raster_dataset.GetGeoTransform()
proj = raster_dataset.GetProjectionRef()

bands_data = []
for b in range(1, raster_dataset.RasterCount+1):
    band = raster_dataset.GetRasterBand(b)
    bands_data.append(band.ReadAsArray())

Sin embargo, cuando intento iterar sobre las bandas, se consume toda la memoria y no se carga durante horas.

¿Alguna sugerencia sobre otras formas de manejar grandes datos espaciales en python? He utilizado ArcGIS y R en el pasado, pero soy relativamente nuevo en python. Preferiría trabajar desde la línea de comandos, si es posible.

Answer 1

Una gran característica de los datos ráster es que suelen permitir el procesamiento en bloques. Puedes "dividir" el ráster en ventanas rectangulares para reducir la huella de memoria de tu proceso, o para procesar bloques en paralelo y obtener resultados más rápidamente.

La documentación de los enlaces Python de GDAL es escasa y los ejemplos de obtención de ventanas rasterizadas mediante ReadAsArray() son escasos, pero he encontrado esto en las pruebas de GDAL: https://github.com/OSGeo/gdal/blob/77544764f51420e42468641fc3d5a087f8ea6d8f/autotest/gcore/numpy_rw.py#L115 . Pasa algunos desplazamientos y la anchura y altura de tu ventana y tendrás un subconjunto de bandas como un array numpy.

Para aplicar esto en su programa, haga un bucle sobre los bloques de la trama, y dentro de este bucle sobre las bandas de la trama.

Answer 2

Cuando ReadAsArray() estás creando un array numpy, lo que significa esencialmente que lo estás cargando en memoria. Para un TIF de 30GB, espera algo más o menos así. Es decir, sólo podrás trabajar con él si tienes una estación de trabajo o un servidor de más de 64 GB. E incluso entonces, no es deseable.

Para rastreos más grandes, considere trabajar con bloques. Aquí hay un buen tutorial sobre esto (a partir de la página 24). Tenga en cuenta, sin embargo, que este documento fue escrito para una versión anterior de GDAL. Desde la versión 2.0, la obtención del tamaño de bloque personalizado de un raster se realiza con el método GDALRasterBand::GetBlockSize().

También hay que tener en cuenta que una trama se itera de forma más eficiente sobre su tamaño de bloque por defecto (cada formato de trama tendrá el suyo propio). Sin embargo, el tipo de bloque de un TIFF puede ser un mosaico o una tira. Mosaico significa un bloque rectangular, y los TIFFs creados con GDAL tienen ese tamaño por defecto (un mosaico de 256x256, específicamente), pero los TIFFs con tipo de tira tienen bloques que abarcan todo el ancho del raster, e incluso pueden tener sólo una fila de altura. Esto no sólo dificulta ciertos procesos (como el kernel), sino que también puede requerir demasiada memoria. Compruebe de antemano el tamaño de los bloques y, si se trata de una tira, elija un tamaño de bloque personalizado.