Инструменты растровой обработки

В этом примере показаны некоторые инструменты обработки растра и изображений, созданные KDRM для автоматизации рабочих процессов с растрами. Эти инструменты полезны для таких задач, как пакетное вырезание растров в интересующую область, повторная выборка растров до требуемого выходного разрешения, подготовка входных растровых данных для последующего анализа и удаление медленных повторяющихся шагов ручной обработки из более крупных рабочих процессов доставки.

Приведенные ниже примеры относятся только к пакетам Python с открытым исходным кодом. На практике подобные инструменты можно запускать как автономные сценарии, запланированные задания обработки или встраивать в более широкие рабочие процессы рабочего стола в ArcGIS Pro или QGIS, в зависимости от клиентской среды и способа управления рабочим процессом.

Для чего используются эти инструменты

• Пакетное вырезание DEM, DTM и других растров в интересующую область.
• Повторная выборка растров для выравнивания размеров ячеек для последующих расчетов и моделирования.
• Фильтрация и проверка входных растровых данных перед последующим анализом.
• Автоматизация этапов обработки растра, которые в противном случае потребовали бы многих часов ручной работы.
• Создание повторяемых конвейеров обработки растра, которые можно запускать в настольных, серверных средах или средах командной строки.

Почему мы создаем подобные инструменты

Мы регулярно создаем собственную систему автоматизации для проектов с большим количеством растровых данных, поскольку производительность, повторяемость и переносимость имеют значение. Иногда требуется запустить рабочий процесс с открытым исходным кодом без зависимости от ArcGIS. В других случаях та же логика должна находиться внутри рабочего процесса ArcGIS Pro или QGIS, чтобы аналитики могли запускать процесс в знакомой среде рабочего стола. Базовый подход тот же: автоматизируйте повторяющиеся части, сохраняйте прозрачность логики и упрощайте повторный запуск и поддержку рабочего процесса.

Приведенный ниже пример обрезки был создан для рабочего процесса, в котором пакетная обрезка растров с использованием традиционного настольного подхода была слишком медленной. Переписывание процесса с помощью пакетов Python с открытым исходным кодом сделало рабочий процесс значительно более быстрым и практичным для больших коллекций растров. Пример повторной выборки был создан по той же причине: большие наборы растровых данных необходимо было выровнять по более мелкому размеру выходной ячейки, чтобы последующие растровые вычисления могли выполняться с требуемым разрешением.

Пример 1. Пакетное вырезание растров с помощью Python с открытым исходным кодом

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

Пример 2. Повторная выборка растров с помощью Python с открытым исходным кодом

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

Варианты доставки

Приведенные выше примеры намеренно просты, но тот же подход можно распространить на более крупные производственные рабочие процессы с ведением журнала, проверкой, параллельной обработкой, проверкой метаданных, упаковкой и интеграцией в более широкие конвейеры доставки ГИС. В зависимости от требований клиента мы можем предоставить эти инструменты в виде автономных сценариев, инструментов, интегрированных в рабочий стол для ArcGIS Pro или QGIS, или готовых сервисов как часть более широкой платформы пространственных данных.