Ferramentas de processamento raster

Este exemplo mostra algumas das ferramentas raster e de processamento de imagem que o KDRM criou para automatizar fluxos de trabalho raster. Essas ferramentas são úteis para tarefas como recorte de lote de rasters para uma área de interesse, reamostragem de rasters para uma resolução de saída necessária, preparação de entradas de raster para análise posterior e remoção de etapas lentas e repetitivas de processamento manual de fluxos de trabalho de entrega maiores.

Os exemplos abaixo concentram-se apenas em pacotes Python de código aberto. Na prática, ferramentas como essa podem ser executadas como scripts independentes, trabalhos de processamento agendados ou incorporadas em fluxos de trabalho de desktop mais amplos no ArcGIS Pro ou QGIS, dependendo do ambiente do cliente e de como o fluxo de trabalho precisa ser operado.

Para que essas ferramentas são usadas

• Recorte em lote DEMs, DTMs e outros rasters para uma área de interesse.
• Reamostragem de rasters para alinhar tamanhos de células para cálculos e modelagem posteriores.
• Filtrar e validar entradas raster antes da análise downstream.
• Automatização de etapas de processamento raster que, de outra forma, levariam horas de trabalho manual.
• Construir pipelines de processamento raster repetíveis que podem ser executados em ambientes de desktop, servidor ou linha de comando.

Por que construímos ferramentas como esta

Desenvolvemos regularmente automação personalizada para projetos com muitos rasters porque o desempenho, a repetibilidade e a portabilidade são importantes. Às vezes, o requisito é executar como um fluxo de trabalho de código aberto puro, sem dependência do ArcGIS. Em outros casos, a mesma lógica precisa estar dentro de um fluxo de trabalho do ArcGIS Pro ou QGIS para que os analistas possam executar o processo em um ambiente de área de trabalho familiar. A abordagem subjacente é a mesma: automatizar as partes repetitivas, manter a lógica transparente e tornar o fluxo de trabalho mais fácil de executar e suportar.

O exemplo de recorte abaixo foi criado para um fluxo de trabalho em que o recorte de rasters em lote com uma abordagem de desktop tradicional era muito lento. Reescrever o processo com pacotes Python de código aberto tornou o fluxo de trabalho muito mais rápido e prático para coleções raster maiores. O exemplo de reamostragem foi construído por um motivo semelhante, onde grandes conjuntos de dados raster precisavam ser alinhados a um tamanho de célula de saída mais preciso para que cálculos raster posteriores pudessem ser executados na resolução necessária.

Exemplo 1. Rasters de recorte em lote com Python de código aberto

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

Exemplo 2. Reamostragem de rasters com Python de código aberto

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

Opções de entrega

Os exemplos acima são deliberadamente simples, mas a mesma abordagem pode ser estendida para fluxos de trabalho de produção maiores com registro, validação, processamento paralelo, verificações de metadados, empacotamento e integração em pipelines de entrega de GIS mais amplos. Dependendo da necessidade do cliente, podemos fornecer essas ferramentas como scripts independentes, ferramentas integradas ao desktop para ArcGIS Pro ou QGIS, ou serviços produzidos como parte de uma plataforma de dados espaciais mais ampla.