أدوات معالجة البيانات النقطية

يوضح هذا المثال بعض أدوات معالجة الصور والنقطية التي أنشأها KDRM لأتمتة سير عمل البيانات النقطية. تعتبر هذه الأدوات مفيدة لمهام مثل قطع البيانات النقطية دفعة واحدة إلى منطقة الاهتمام، وإعادة تشكيل البيانات النقطية إلى دقة الإخراج المطلوبة، وإعداد مدخلات البيانات النقطية للتحليل اللاحق، وإزالة خطوات المعالجة اليدوية البطيئة والمتكررة من سير عمل التسليم الأكبر.

تركز الأمثلة أدناه على حزم Python مفتوحة المصدر فقط. من الناحية العملية، يمكن تشغيل أدوات مثل هذه كبرامج نصية مستقلة أو مهام معالجة مجدولة أو مضمنة داخل مسارات عمل سطح المكتب الأوسع في ArcGIS Pro أو QGIS، اعتمادًا على بيئة العميل وكيفية تشغيل سير العمل.

فيما تستخدم هذه الأدوات

• قطع دفعة من نماذج DEM وDTMs وغيرها من البيانات النقطية إلى منطقة الاهتمام.
• إعادة تشكيل البيانات النقطية لمحاذاة أحجام الخلايا لإجراء العمليات الحسابية والنمذجة اللاحقة.
• تصفية المدخلات النقطية والتحقق من صحتها قبل التحليل النهائي.
• أتمتة خطوات معالجة البيانات النقطية التي قد تستغرق ساعات من العمل اليدوي.
• إنشاء خطوط أنابيب معالجة نقطية قابلة للتكرار يمكن تشغيلها في بيئات سطح المكتب أو الخادم أو سطر الأوامر.

لماذا نبني أدوات مثل هذه

نحن نبني بشكل منتظم أتمتة مخصصة للمشاريع ذات البيانات النقطية الثقيلة نظرًا لأهمية الأداء وقابلية التكرار وقابلية النقل. في بعض الأحيان يكون المتطلب هو التشغيل كسير عمل خالص مفتوح المصدر بدون تبعية 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، أو خدمات إنتاجية كجزء من نظام أساسي أوسع للبيانات المكانية.