¿Cómo fusionar rásters con rasterio en bloques para evitar MemoryError?

Question

Estoy tratando de fusionar dos rasters con rasterio.

Sin embargo, cuando ejecuto dest_array, out_transform = merge(sources), recibo el siguiente error:

MemoryError: No se puede asignar matriz con forma (1, 56160, 85196) y tipo de datos float64

Sé que rasterio puede leer y escribir archivos en bloques.

Desafortunadamente, no parece haber una forma de hacerlo al llamar a merge.

¿Hay alguna manera de usar bloques con merge u otra forma de fusionar rasters sin usar el merge método?

Sé cómo leer y escribir datos en bloques, pero no tengo idea de cómo hacerlo al usar el método merge.

Aquí está el código que estoy usando hasta la línea que llama a merge y que me da el error:

import rasterio
from rasterio.merge import merge
import glob
import os

def merge_raster_in_folder(src_dir, dest_dir):
    """Fusionar todos los rasters (extensión TIFF) dentro del directorio dado"""

    # crear lista de archivos de raster en la carpeta de origen
    rasters = []
    os.chdir(src_dir)
    for file in glob.glob("*.tif"):
        rasters.append(file)

    # crear lista de objetos de raster a partir de la lista de nombres de archivo de raster
    sources = [rasterio.open(raster) for raster in rasters]

    # crear la matriz que representa todos los rasters de origen unidos
    dest_array, out_transform = merge(sources)

Answer 1

Finalmente resolví el problema usando numpy.memmap para

crear un mapa de memoria a un array almacenado en un archivo binario en disco

y luego procesando los rásteres de entrada en ventanas y bloques. Puede ser más lento pero funciona y estoy contento con el resultado (necesito agradecer al usuario @Thomas que me ayudó en algunos pasos).

El código que estoy usando está tomado y modificado del código fuente de rasterio merge.py, hasta la parte donde crea un "enorme" (en mi caso) array de numpy de ceros que estaba causando un MemoryError en mi caso.

Aquí está el código final:

from tempfile import mkdtemp
import rasterio
from rasterio import Affine
from rasterio import windows
import math
import numpy as np
import os

INPUT_FILES = [r'ruta/a/raster1.tif', r'ruta/a/rasterN.tif']

sources = [rasterio.open(raster) for raster in INPUT_FILES]

memmap_file = os.path.join(mkdtemp(), 'test.mymemmap')

# adaptado de https://github.com/mapbox/rasterio/blob/master/rasterio/merge.py
first = sources[0]
first_res = first.res
dtype = first.dtypes[0]
# Determinar cantidad de bandas de salida
output_count = first.count

# Extensión de todas las entradas
# escanear archivos de entrada
xs = []
ys = []
for src in sources:
    left, bottom, right, top = src.bounds
    xs.extend([left, right])
    ys.extend([bottom, top])
dst_w, dst_s, dst_e, dst_n = min(xs), min(ys), max(xs), max(ys)

out_transform = Affine.translation(dst_w, dst_n)

# Resolución/tamaño de píxel
res = first_res
out_transform *= Affine.scale(res[0], -res[1])

# Calcular forma del array de salida. Garantizamos que cubrirá completamente los límites de salida
output_width = int(math.ceil((dst_e - dst_w) / res[0]))
output_height = int(math.ceil((dst_n - dst_s) / res[1]))

# Ajustar límites para que encajen
dst_e, dst_s = out_transform * (output_width, output_height)

# crear array de destino
# forma del array de destino
shape = (output_height, output_width)
# dest = np.zeros((output_count, output_height, output_width), dtype=dtype)
# Usando numpy.memmap para crear arrays directamente mapeados en un archivo
dest_array = np.memmap(memmap_file, dtype=dtype,
                       mode='w+', shape=shape)

dest_profile = {
    "driver": 'GTiff',
    "height": dest_array.shape[0],
    "width": dest_array.shape[1],
    "count": output_count,
    "dtype": dest_array.dtype,
    "crs": '+proj=latlong',
    "transform": out_transform
}

# abrir archivo de salida en modo escritura/lectura y completar con el array de mosaico de destino
with rasterio.open(
    os.path.join(r'out/dir', 'test.tif'),
    'w+',
    **dest_profile
) as mosaic_raster:
    for src in sources:
        for ji, src_window in src.block_windows(1):

            print(ji)
            r = src.read(1, window=src_window)
            # almacenar valor de nodata del ráster
            nodata = src.nodatavals[0]
            # reemplazar ceros con nan
            r[r == nodata] = np.nan
            # convertir ubicación de ventana de entrada relativa a ubicación de ventana de salida relativa
            # usando coordenadas del mundo real (límites)
            src_bounds = windows.bounds(
                src_window, transform=src.profile["transform"])
            dst_window = windows.from_bounds(
                *src_bounds, transform=mosaic_raster.profile["transform"])

            # redondear los valores de dest_window ya que pueden ser decimales
            dst_window = windows.Window(round(dst_window.col_off), round(
                dst_window.row_off), round(dst_window.width), round(dst_window.height))
            # antes de escribir la ventana, reemplazar nodata de origen con nodata de destino ya que podría haber sido escrito previamente (por ejemplo, otro país adyacente)
            # https://stackoverflow.com/a/43590909/1979665
            dest_pre = mosaic_raster.read(1, window=dst_window)
            mask = (np.isnan(r))
            r_mod = np.copy(r)
            r_mod[mask] = dest_pre[mask]
            mosaic_raster.write(r_mod, 1, window=dst_window)

os.remove(memmap_file)

Answer 2

Me enfrenté al mismo problema. La solución más fácil, en este caso, es no usar rasterio. En su lugar, usar gdal será mucho más fácil, sin necesidad de modificar el código fuente. Siga el segundo paso de este video cuidadosamente.

Básicamente, están creando un archivo XML de ráster virtual utilizando la función gdal.BuildVRT() de gdal. El archivo XML de ráster virtual almacena todos los rásteres individuales que se van a fusionar. Y la función gdal.Translate() se utiliza para convertir el XML a un GTiff. La función BuildVRT() evita el uso de la función numpy.zeros() en rasterio.merge(), que es la fuente de error de memoria en el caso de conjuntos de datos grandes.