Alat Pemrosesan Raster

Contoh ini menunjukkan beberapa alat raster dan pemrosesan gambar yang dibuat KDRM untuk mengotomatisasi alur kerja raster. Alat-alat ini berguna untuk tugas-tugas seperti memotong raster secara batch ke area yang diminati, mengambil sampel ulang raster ke resolusi keluaran yang diperlukan, menyiapkan masukan raster untuk analisis selanjutnya, dan menghilangkan langkah-langkah pemrosesan manual berulang yang lambat dari alur kerja pengiriman yang lebih besar.

Contoh di bawah ini hanya fokus pada paket open-source Python saja. Dalam praktiknya, alat seperti ini dapat dijalankan sebagai skrip mandiri, tugas pemrosesan terjadwal, atau tertanam dalam alur kerja desktop yang lebih luas di ArcGIS Pro atau QGIS, bergantung pada lingkungan klien dan cara alur kerja tersebut dioperasikan.

Untuk apa alat ini digunakan

• Batch kliping DEM, DTM, dan raster lainnya ke area yang diinginkan.
• Pengambilan sampel ulang raster untuk menyelaraskan ukuran sel untuk penghitungan dan pemodelan selanjutnya.
• Memfilter dan memvalidasi input raster sebelum analisis hilir.
• Mengotomatiskan langkah-langkah pemrosesan raster yang membutuhkan waktu berjam-jam kerja manual.
• Membangun jalur pemrosesan raster berulang yang dapat dijalankan di lingkungan desktop, server, atau baris perintah.

Mengapa kami membuat alat seperti ini

Kami secara teratur membangun otomatisasi khusus untuk proyek-proyek raster-berat karena kinerja, pengulangan, dan portabilitas penting. Terkadang persyaratannya adalah dijalankan sebagai alur kerja sumber terbuka murni tanpa ketergantungan ArcGIS. Dalam kasus lain, logika yang sama perlu ditempatkan di dalam alur kerja ArcGIS Pro atau QGIS sehingga analis dapat menjalankan proses dalam lingkungan desktop yang familiar. Pendekatan dasarnya sama: mengotomatiskan bagian yang berulang, menjaga logika tetap transparan, dan membuat alur kerja lebih mudah untuk dijalankan kembali dan didukung.

Contoh kliping di bawah ini dibuat untuk alur kerja di mana kliping batch raster dengan pendekatan desktop tradisional terlalu lambat. Menulis ulang proses dengan paket Python sumber terbuka membuat alur kerja jauh lebih cepat dan praktis untuk koleksi raster yang lebih besar. Contoh pengambilan sampel ulang dibuat untuk alasan serupa, yaitu kumpulan data raster yang besar perlu disejajarkan dengan ukuran sel keluaran yang lebih halus sehingga penghitungan raster selanjutnya dapat berjalan pada resolusi yang diperlukan.

Contoh 1. Memotong raster secara batch dengan Python sumber terbuka

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

Contoh 2. Pengambilan sampel ulang raster dengan Python sumber terbuka

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

Opsi pengiriman

Contoh di atas sengaja dibuat sederhana, namun pendekatan yang sama dapat diperluas untuk alur kerja produksi yang lebih besar dengan logging, validasi, pemrosesan paralel, pemeriksaan metadata, pengemasan, dan integrasi ke dalam jalur pengiriman GIS yang lebih luas. Bergantung pada kebutuhan klien, kami dapat menyediakan alat-alat ini sebagai skrip mandiri, alat terintegrasi desktop untuk ArcGIS Pro atau QGIS, atau layanan produksi sebagai bagian dari platform data spasial yang lebih luas.