Herramientas de procesamiento raster

Este ejemplo muestra algunas de las herramientas de procesamiento raster e imagen que KDRM ha creado para automatizar flujos de trabajo raster. Estas herramientas son útiles para tareas como recortar rasters en lote a un área de interés, remuestrear rasters a una resolución de salida requerida, preparar entradas raster para análisis posteriores y eliminar pasos manuales lentos y repetitivos de flujos de trabajo de mayor escala.

Los ejemplos siguientes se centran únicamente en paquetes Python de código abierto. En la práctica, herramientas como estas pueden ejecutarse como scripts independientes, trabajos programados o integrarse dentro de flujos de escritorio en ArcGIS Pro o QGIS, según el entorno del cliente y la forma en que se necesite operar el proceso.

Para qué se utilizan estas herramientas

• Recorte en lote de DEM, DTM y otros rasters a un área de interés.
• Remuestreo de rasters para alinear tamaños de celda antes de cálculos y modelado posteriores.
• Filtrado y validación de entradas raster antes del análisis downstream.
• Automatización de pasos raster que de otro modo requerirían horas de trabajo manual.
• Construcción de canalizaciones raster repetibles que puedan ejecutarse en entornos de escritorio, servidor o línea de comandos.

Por qué construimos herramientas como estas

Construimos automatización personalizada para proyectos con mucho raster porque el rendimiento, la repetibilidad y la portabilidad importan. A veces el requisito es ejecutar un flujo puramente de código abierto sin dependencia de ArcGIS. En otros casos, la misma lógica debe integrarse en un flujo de ArcGIS Pro o QGIS para que los analistas puedan ejecutar el proceso dentro de un entorno de escritorio familiar. El enfoque subyacente es el mismo: automatizar lo repetitivo, mantener la lógica transparente y facilitar que el flujo pueda repetirse y mantenerse.

El ejemplo de recorte de abajo se creó para un flujo donde el recorte masivo de rasters con un enfoque tradicional de escritorio era demasiado lento. Reescribir el proceso con paquetes Python de código abierto hizo el flujo mucho más rápido y práctico para colecciones raster grandes. El ejemplo de remuestreo se construyó por una razón similar, cuando varios datasets raster necesitaban alinearse a un tamaño de celda de salida más fino para que los cálculos posteriores pudieran ejecutarse con la resolución requerida.

Ejemplo 1. Recorte de rasters en lote con Python de código abierto

import json
import numpy as np
import os
import rasterio
import shutil

from rasterio.mask import mask
from shapely.geometry import shape
# the input and output directories for the rasters
input_directory = "input"
output_directory = "output"

# delete and recreate the output directory
if(os.path.exists(output_directory)):
    shutil.rmtree(output_directory)

os.makedirs(output_directory)

feature_collection = dict()

with open('national_park.geojson', 'r') as f:
    feature_collection = json.load(f)
area_of_interest = feature_collection['features'][0]['geometry']
clip_area = shape(area_of_interest)

for file in os.listdir(input_directory):
    # filter by tif files and only files with 2019 in the name
    if file.endswith('.tif') and '2019' in file:
        print(file)

        input_raster = os.path.join(input_directory, file)
        output_raster=os.path.join(output_directory, file)
        with rasterio.open(input_raster) as src:
            try:
                out_image, out_transform = mask(src, [clip_area], crop=True)

                # Set all nodata values to -99999
                out_image[out_image == src.nodata] = -99999
                raster_intersected = not np.all(out_image == -99999)
                # Check if the clip area intersects the raster and only save the raster if the clip area intersects with it
                if raster_intersected:
                    out_meta = src.meta.copy()
                    out_meta.update(
                        {
                            "driver": "GTiff",
                            "height": out_image.shape[1],
                            "width": out_image.shape[2],
                            "transform": out_transform,
                        }
                    )
                    with rasterio.open(output_raster, "w", **out_meta) as dest:
                        dest.write(out_image)

            except Exception as err:
                if ('Input shapes do not overlap raster' not in str(err)):
                    print('### Error with {} ###'.format(input_raster))
                    print(err)
                pass

Ejemplo 2. Remuestreo de rasters con Python de código abierto

import os
from osgeo import gdal

import numpy as np
import os

# the input and output directories
input_directory = 'input/DTM'
output_directory = 'output/DTM'

if not os.path.exists(output_directory):
    os.makedirs(output_directory)
for root, dirs, files in os.walk(input_directory):
    for file in files:
        if(file.endswith('.tif')):
            input_file = os.path.join(root, file)
            output_file = os.path.join(output_directory, file)

            print(file)

            # Read the input raster file
            raster = gdal.Open(input_file)
            # Get the original cell size
            cell_size = raster.GetGeoTransform()[1]
            # print('Original cell size: {}'.format(cell_size))

            # Calculate the number of cells for the output
            num_cells = int(np.ceil(cell_size / 0.5))
            # Create the output raster file
            driver = gdal.GetDriverByName('GTiff')
            out_raster = driver.Create(output_file, raster.RasterXSize * num_cells,
                                    raster.RasterYSize * num_cells, 1,
                                    gdal.GDT_Float32)

            # Set the output raster's projection
            out_raster.SetProjection(raster.GetProjection())
            # Set the output raster's geotransform
            out_raster_geo = list(raster.GetGeoTransform())
            out_raster_geo[1] = 0.5
            out_raster_geo[5] = -0.5
            out_raster.SetGeoTransform(out_raster_geo)

            nodata = -99999

            for i in range(1, out_raster.RasterCount + 1):
                # set the nodata value of the band
                out_raster.GetRasterBand(i).SetNoDataValue(nodata)
            # Copy the input raster data to the output raster
            gdal.ReprojectImage(raster, out_raster,
                                raster.GetProjection(),
                                out_raster.GetProjection(),
                                gdal.GRA_NearestNeighbour)

            # Close the files
            out_raster = None
            raster = None

Opciones de entrega

Los ejemplos anteriores son deliberadamente simples, pero el mismo enfoque puede ampliarse para flujos de producción más grandes con logging, validación, procesamiento en paralelo, controles de metadatos, empaquetado e integración en canalizaciones GIS más amplias. Según el requisito del cliente, podemos proporcionar estas herramientas como scripts independientes, herramientas integradas en ArcGIS Pro o QGIS, o servicios productivizados dentro de una plataforma espacial más amplia.