Strumenti di elaborazione raster

Questo esempio mostra alcuni degli strumenti raster e di elaborazione delle immagini creati da KDRM per automatizzare i flussi di lavoro raster. Questi strumenti sono utili per attività quali il ritaglio batch di raster in un’area di interesse, il ricampionamento dei raster alla risoluzione di output richiesta, la preparazione di input raster per un’analisi successiva e la rimozione di passaggi lenti e ripetitivi di elaborazione manuale da flussi di lavoro di distribuzione più ampi.

Gli esempi seguenti si concentrano solo sui pacchetti Python open source. In pratica, strumenti come questo possono essere eseguiti come script autonomi, processi di elaborazione pianificati o incorporati all’interno di flussi di lavoro desktop più ampi in ArcGIS Pro o QGIS, a seconda dell’ambiente client e di come deve essere gestito il flusso di lavoro.

A cosa servono questi strumenti

• Ritaglio batch di DEM, DTM e altri raster in un’area di interesse.
• Ricampionamento dei raster per allineare le dimensioni delle celle per calcoli e modellazioni successivi.
• Filtraggio e convalida degli input raster prima dell’analisi downstream.
• Automatizzazione delle fasi di elaborazione raster che altrimenti richiederebbero ore di lavoro manuale.
• Creazione di pipeline di elaborazione raster ripetibili che possono essere eseguite in ambienti desktop, server o a riga di comando.

Perché creiamo strumenti come questo

Realizziamo regolarmente automazioni personalizzate per progetti ad uso intensivo di raster perché prestazioni, ripetibilità e portabilità sono importanti. A volte il requisito è quello di eseguirlo come un flusso di lavoro open source puro senza una dipendenza ArcGIS. In altri casi la stessa logica deve essere inserita in un flusso di lavoro ArcGIS Pro o QGIS in modo che gli analisti possano eseguire il processo in un ambiente desktop familiare. L’approccio di fondo è lo stesso: automatizzare le parti ripetitive, mantenere la logica trasparente e rendere il flusso di lavoro più facile da rieseguire e supportare.

L’esempio di ritaglio riportato di seguito è stato creato per un flusso di lavoro in cui il ritaglio batch di raster con un approccio desktop tradizionale era troppo lento. La riscrittura del processo con pacchetti Python open source ha reso il flusso di lavoro notevolmente più veloce e più pratico per raccolte raster più grandi. L’esempio di ricampionamento è stato creato per un motivo simile, in cui i set di dati raster di grandi dimensioni dovevano essere allineati a una dimensione della cella di output più precisa in modo che i calcoli raster successivi potessero essere eseguiti alla risoluzione richiesta.

Esempio 1. Ritaglio batch di raster con Python open source

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

Esempio 2. Ricampionamento dei raster con Python open source

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

Opzioni di consegna

Gli esempi sopra riportati sono volutamente semplici, ma lo stesso approccio può essere esteso a flussi di lavoro di produzione più ampi con registrazione, convalida, elaborazione parallela, controlli dei metadati, confezionamento e integrazione in pipeline di distribuzione GIS più ampie. A seconda delle esigenze del cliente, possiamo fornire questi strumenti come script autonomi, strumenti integrati nel desktop per ArcGIS Pro o QGIS o servizi produttivi come parte di una piattaforma di dati spaziali più ampia.