Raster İşleme Araçları

Bu örnek, KDRM’nin tarama iş akışlarını otomatikleştirmek için oluşturduğu tarama ve görüntü işleme araçlarından bazılarını göstermektedir. Bu araçlar, rasterleri ilgi alanına toplu olarak kırpmak, rasterleri gerekli çıktı çözünürlüğüne göre yeniden örneklemek, tarama girdilerini daha sonraki analiz için hazırlamak ve daha büyük teslim iş akışlarından yavaş tekrarlanan manuel işleme adımlarını kaldırmak gibi görevler için kullanışlıdır.

Aşağıdaki örnekler yalnızca açık kaynaklı Python paketlerine odaklanmaktadır. Uygulamada bunun gibi araçlar, istemci ortamına ve iş akışının nasıl çalıştırılması gerektiğine bağlı olarak bağımsız komut dosyaları, zamanlanmış işleme işleri olarak çalıştırılabilir veya ArcGIS Pro veya QGIS’teki daha geniş masaüstü iş akışlarına yerleştirilebilir.

Bu araçlar ne için kullanılıyor?

• DEM’leri, DTM’leri ve diğer taramaları ilgilenilen bir alana toplu olarak kırpın.
• Daha sonraki hesaplamalar ve modelleme için hücre boyutlarını hizalamak üzere rasterleri yeniden örnekleme.
• Aşağı akış analizinden önce tarama girdilerinin filtrelenmesi ve doğrulanması.
• Aksi takdirde saatlerce manuel çalışma gerektirecek tarama işleme adımlarının otomatikleştirilmesi.
• Masaüstünde, sunucuda veya komut satırı ortamlarında çalıştırılabilen tekrarlanabilir tarama işleme hatları oluşturmak.

Neden bunun gibi araçlar geliştiriyoruz?

Performans, tekrarlanabilirlik ve taşınabilirlik önemli olduğundan tarama ağırlıklı projeler için düzenli olarak özel otomasyon geliştiriyoruz. Bazen gereksinim, ArcGIS bağımlılığı olmadan saf bir açık kaynak iş akışı olarak çalışmaktır. Diğer durumlarda, analistlerin süreci tanıdık bir masaüstü ortamında yürütebilmesi için aynı mantığın ArcGIS Pro veya QGIS iş akışında da bulunması gerekir. Temel yaklaşım aynıdır: Tekrarlanan parçaları otomatikleştirin, mantığı şeffaf tutun ve iş akışının yeniden çalıştırılmasını ve desteklenmesini kolaylaştırın.

Aşağıdaki kırpma örneği, geleneksel masaüstü yaklaşımıyla taramaların toplu olarak kırpılmasının çok yavaş olduğu bir iş akışı için oluşturulmuştur. Sürecin açık kaynaklı Python paketleriyle yeniden yazılması, iş akışını önemli ölçüde daha hızlı ve daha büyük tarama koleksiyonları için daha pratik hale getirdi. Yeniden örnekleme örneği de benzer bir nedenden ötürü oluşturulmuştur; daha sonraki tarama hesaplamalarının gerekli çözünürlükte çalışabilmesi için büyük tarama veri kümelerinin daha hassas bir çıktı hücresi boyutuna hizalanması gerekiyordu.

Örnek 1. Açık kaynak Python ile rasterleri toplu kırpma

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

Örnek 2. Açık kaynak Python ile rasterleri yeniden örnekleme

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

Teslimat seçenekleri

Yukarıdaki örnekler kasıtlı olarak basittir ancak aynı yaklaşım, günlüğe kaydetme, doğrulama, paralel işleme, meta veri kontrolleri, paketleme ve daha geniş GIS dağıtım hatlarına entegrasyon ile daha büyük üretim iş akışları için genişletilebilir. Müşteri ihtiyacına bağlı olarak bu araçları bağımsız komut dosyaları, ArcGIS Pro veya QGIS için masaüstüne entegre araçlar veya daha geniş bir mekansal veri platformunun parçası olarak üretime geçirilmiş hizmetler olarak sağlayabiliriz.