Công cụ xử lý raster
Ví dụ này cho thấy một số công cụ xử lý hình ảnh và raster mà KDRM đã tạo để tự động hóa quy trình làm việc raster. Những công cụ này rất hữu ích cho các tác vụ như cắt các raster hàng loạt đến một khu vực quan tâm, lấy mẫu lại các raster để có độ phân giải đầu ra cần thiết, chuẩn bị đầu vào raster để phân tích sau này và loại bỏ các bước xử lý thủ công lặp đi lặp lại chậm khỏi quy trình phân phối lớn hơn.
Các ví dụ bên dưới chỉ tập trung vào các gói Python nguồn mở. Trong thực tế, các công cụ như thế này có thể chạy dưới dạng tập lệnh độc lập, công việc xử lý theo lịch trình hoặc được nhúng bên trong quy trình làm việc trên máy tính để bàn rộng hơn trong ArcGIS Pro hoặc QGIS, tùy thuộc vào môi trường máy khách và cách vận hành quy trình công việc.
Những công cụ này dùng để làm gì
- Cắt hàng loạt DEM, DTM và các trình quét khác vào khu vực quan tâm.
- Lấy mẫu lại các raster để căn chỉnh kích thước ô cho các tính toán và mô hình hóa sau này.
- Lọc và xác nhận đầu vào raster trước khi phân tích tiếp theo.
- Tự động hóa các bước xử lý raster mà nếu không sẽ phải mất hàng giờ làm việc thủ công.
- Xây dựng các quy trình xử lý raster lặp lại có thể chạy trong môi trường máy tính để bàn, máy chủ hoặc dòng lệnh.
Tại sao chúng tôi xây dựng những công cụ như thế này
Chúng tôi thường xuyên xây dựng tính năng tự động hóa tùy chỉnh cho các dự án nặng về raster vì hiệu suất, độ lặp lại và tính di động rất quan trọng. Đôi khi, yêu cầu là chạy như một quy trình làm việc nguồn mở thuần túy mà không phụ thuộc vào ArcGIS. Trong các trường hợp khác, logic tương tự cần phải có trong quy trình làm việc của ArcGIS Pro hoặc QGIS để các nhà phân tích có thể chạy quy trình trong môi trường máy tính để bàn quen thuộc. Cách tiếp cận cơ bản là giống nhau: tự động hóa các phần lặp đi lặp lại, giữ logic minh bạch và làm cho quy trình làm việc dễ dàng chạy lại và hỗ trợ hơn.
Ví dụ cắt bớt bên dưới được tạo cho một quy trình làm việc trong đó việc cắt các trình quét hàng loạt bằng cách tiếp cận trên máy tính để bàn truyền thống quá chậm. Viết lại quy trình bằng các gói Python nguồn mở giúp quy trình làm việc nhanh hơn đáng kể và thiết thực hơn cho các bộ sưu tập raster lớn hơn. Ví dụ về việc lấy mẫu lại được xây dựng vì một lý do tương tự, trong đó các bộ dữ liệu raster lớn cần được căn chỉnh theo kích thước ô đầu ra mịn hơn để các phép tính raster sau này có thể chạy ở độ phân giải cần thiết.
Ví dụ 1. Cắt hàng loạt các trình quét bằng Python nguồn mở
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)
passVí dụ 2. Lấy mẫu lại các trình quét bằng Python nguồn mở
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 = NoneTùy chọn giao hàng
Các ví dụ trên tuy đơn giản nhưng có thể mở rộng cách tiếp cận tương tự cho các quy trình sản xuất lớn hơn với tính năng ghi nhật ký, xác thực, xử lý song song, kiểm tra siêu dữ liệu, đóng gói và tích hợp vào các quy trình phân phối GIS rộng hơn. Tùy thuộc vào yêu cầu của khách hàng, chúng tôi có thể cung cấp các công cụ này dưới dạng tập lệnh độc lập, công cụ tích hợp trên máy tính để bàn cho ArcGIS Pro hoặc QGIS hoặc các dịch vụ được sản xuất như một phần của nền tảng dữ liệu không gian rộng hơn.