光柵處理工具

此範例展示了 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

交付選項

上面的範例故意簡單，但相同的方法可以擴展到更大的生產工作流程，包括日誌記錄、驗證、並行處理、元資料檢查、打包以及整合到更廣泛的 GIS 交付管道中。根據客戶要求，我們可以將這些工具作為獨立腳本、ArcGIS Pro 或 QGIS 桌面整合工具或作為更廣泛空間資料平台一部分的生產化服務提供。