Hướng dẫn geocomputation with python - tính toán địa lý với python

Giới thiệu

Trang web này chứa các ý tưởng, mã và một phác thảo của một cuốn sách chưa từng viết về tính toán địa lý với Python.

Động lực

Geocompting với Python, được thúc đẩy bởi sự cần thiết của một nguồn tài nguyên giới thiệu nhưng nghiêm ngặt và được duy trì khi làm việc với dữ liệu địa lý trong Python. Một tính năng độc đáo của cuốn sách là nó thể hiện mã để làm việc với cả các loại dữ liệu địa lý vector và raster. Có nhiều tài nguyên trên các gói Python cho nghiên cứu địa lý và các ứng dụng khác nhau, nhưng theo hiểu biết tốt nhất của chúng tôi, không có tài nguyên nào khác tập hợp các tính năng sau vào một ngôi nhà duy nhất:

  1. Sách giáo khoa giới thiệu nhỏ tập trung vào việc thực hiện tốt các hoạt động cơ bản
  2. Tích hợp các bộ dữ liệu vector và raster trong cùng một cuốn sách và trong mỗi phần
  3. Giải thích rõ ràng về mã và bài tập để tối đa hóa việc học cho người mới đến
  4. Cung cấp các bộ dữ liệu ví dụ sáng suốt và các hoạt động có ý nghĩa để minh họa bản chất ứng dụng của nghiên cứu địa lý

Cuốn sách nhằm mục đích bổ sung các tài nguyên khác trong hệ sinh thái, như được nhấn mạnh bằng cách so sánh với phạm vi cuốn sách với các tác phẩm hiện có và đang tiến hành:

  • Học phân tích không gian địa lý với Python và xử lý địa lý với Python tập trung vào việc xử lý dữ liệu không gian bằng giao diện Python cấp thấp cho GDAL, chẳng hạn như các gói
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    0,
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    1 và
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    2 từ
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    3. Cách tiếp cận này phức tạp hơn, ít pythonic, và có lẽ đã lỗi thời trong việc phát triển các gói như
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    4 và
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    5 được đề cập ở đây
  • Pythongis.org (ở giai đoạn phát triển đầu tiên) tìm cách giới thiệu chung về ‘GIS trong Python, với các phần tập trung vào các yếu tố cần thiết của Python, sử dụng Python với GIS và nghiên cứu trường hợp. So với pythongis.org, GeoCompy có phạm vi tương đối hẹp (1) và tập trung nhiều hơn vào khả năng tương tác vector raster
  • GeographicData.Science là một dự án đầy tham vọng với các chương dành riêng cho các chủ đề nâng cao, với Chương 4 về trọng số không gian đi vào các chủ đề phức tạp tương đối sớm, chẳng hạn. GeoCompy sẽ ngắn hơn, đơn giản hơn và giới thiệu hơn, và bao gồm dữ liệu raster và vector với tầm quan trọng như nhau (1 đến 4)

GeoCompy là một dự án chị em về tính toán địa lý với R - một cuốn sách về phân tích dữ liệu địa lý, trực quan hóa và mô hình hóa bằng ngôn ngữ lập trình R.

Tái tạo cuốn sách này

Một khía cạnh quan trọng của nghiên cứu khoa học và ‘khoa học công dân, đó là sự tham gia là khả năng tái tạo kết quả.

Chất kết dính

Để tái tạo cuốn sách này, bạn chỉ cần nhấp vào liên kết bên dưới để xem mã đang chạy trong trình duyệt web của bạn (xem chi tiết về cách thức hoạt động tại mybinder.org):

Hướng dẫn geocomputation with python - tính toán địa lý với python

Địa phương với quarto

Để chạy mã cục bộ, được khuyến nghị sử dụng tài liệu trên dữ liệu thực, bạn cần có một máy tính hợp lý, ví dụ: & nbsp; với RAM 8 GB. Bạn sẽ cần quyền hành chính để cài đặt các yêu cầu, bao gồm:

  • Một môi trường phát triển tích hợp phù hợp (IDE) như mã vs, rstudio hoặc jupyter máy tính xách tay
  • Quarto, nếu bạn muốn tái tạo trang web truy cập mở sách
  • Một môi trường giống như Anaconda (chúng tôi khuyên bạn nên
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    6) hoặc Docker để có được sự phụ thuộc của hệ thống

Xem dự án READ README để biết chi tiết về việc thiết lập. Sau khi bạn đã cài đặt các phụ thuộc cần thiết và nhân bản hoặc giải nén mã nguồn sách, bạn sẽ có thể sao chép toàn bộ mã với lệnh sau:

Nếu bạn thấy đầu ra như thế bên dưới (với IDE và trình duyệt được sắp xếp để xem các bản cập nhật trực tiếp sau khi chỉnh sửa mã nguồn), xin chúc mừng, nó đã hoạt động!

Hướng dẫn geocomputation with python - tính toán địa lý với python

Địa phương với Jupyter

Ngoài ra, bạn có thể tải xuống và giải nén mã nguồn sách. Thư mục giải nén

from shapely.geometry import Point, LineString, Polygon
points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
line = LineString([(0.4,0.2), (1,0.5)])
poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])

p = gpd.GeoSeries(poly).plot(color='yellow')
p = gpd.GeoSeries(line).plot(ax=p, color='red')
points.plot(ax=p)
7 chứa:

  • Các tệp
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    8 nguồn, một cho mỗi chương
  • Phân hướng phụ
    from shapely.geometry import Point, LineString, Polygon
    points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
    line = LineString([(0.4,0.2), (1,0.5)])
    poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])
    
    p = gpd.GeoSeries(poly).plot(color='yellow')
    p = gpd.GeoSeries(line).plot(ax=p, color='red')
    points.plot(ax=p)
    9 với dữ liệu mẫu được sử dụng trong các phần mã

Giả sử rằng tất cả các gói yêu cầu được cài đặt (xem đầu của mỗi chương), bạn có thể thực thi các tệp

from shapely.geometry import Point, LineString, Polygon
points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
line = LineString([(0.4,0.2), (1,0.5)])
poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])

p = gpd.GeoSeries(poly).plot(color='yellow')
p = gpd.GeoSeries(line).plot(ax=p, color='red')
points.plot(ax=p)
8 thông qua môi trường làm việc đã chọn của bạn (VSCode, Jupyter Notebook, v.v.).

Bài tập

Viết một hàm chấp nhận và mảng và

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
78 chỉ định số lượng hàng/cột để xóa dọc theo một mảng. Hàm cần trả về mảng sửa đổi với các giá trị
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
23 dọc theo các cạnh của nó.

import geopandas as gpd
import numpy as np
import os
import rasterio
import scipy.ndimage

from rasterio.plot import show

Điều kiện tiên quyết

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')

Gói…

Hãy để chúng tôi tải dữ liệu mẫu cho chương này:

Giới thiệu

Hoạt động không gian trên dữ liệu vectorSection 3.3.1), subsets of

0     True
1    False
2     True
dtype: bool
1s can be created with square bracket (
0     True
1    False
2     True
dtype: bool
2) operator using the syntax
0     True
1    False
2     True
dtype: bool
3, where
0     True
1    False
2     True
dtype: bool
4 is an
0     True
1    False
2     True
dtype: bool
1 from which a subset of rows/features will be returned, and
0     True
1    False
2     True
dtype: bool
6 is the ‘subsetting object’.
0     True
1    False
2     True
dtype: bool
6, in turn, can be created using one of the binary geometry relation methods, such as
0     True
1    False
2     True
dtype: bool
8 (see Section 4.3.2).

Để chứng minh tập hợp không gian, chúng tôi sẽ sử dụng các lớp

0     True
1    False
2     True
dtype: bool
9 và
0    False
1    False
2     True
dtype: bool
0, chứa dữ liệu địa lý trên 16 khu vực chính và 101 điểm cao nhất ở New Zealand, tương ứng (Hình & NBSP; 4.1), trong một hệ tọa độ dự kiến. Mã mã sau đây tạo ra một đối tượng đại diện cho Canterbury, sau đó sử dụng tập hợp không gian để trả về tất cả các điểm cao trong khu vực:4.1), in a projected coordinate system. The following code chunk creates an object representing Canterbury, then uses spatial subsetting to return all high points in the region:

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");

Hướng dẫn geocomputation with python - tính toán địa lý với python

Tập hợp không gian 2

non_canterbury_height = nz_height.overlay(canterbury, how='difference')

Kịch bản…

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");

Hướng dẫn geocomputation with python - tính toán địa lý với python

Quan hệ tôpô

from shapely.geometry import Point, LineString, Polygon
points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
line = LineString([(0.4,0.2), (1,0.5)])
poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])

p = gpd.GeoSeries(poly).plot(color='yellow')
p = gpd.GeoSeries(line).plot(ax=p, color='red')
points.plot(ax=p)

Hướng dẫn geocomputation with python - tính toán địa lý với python

0     True
1    False
2     True
dtype: bool

0    False
1    False
2     True
dtype: bool

0     True
1    False
2    False
dtype: bool

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
0

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
1

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
2

Quan hệ tôpô

Quan hệ tôpô

DE-9IM chuỗiSection 3.3.3. Spatial data joining applies the same concept, but instead relies on spatial relations, described in the previous section. As with attribute data, joining adds new columns to the target object (the argument x in joining functions), from a source object (y).

Tham gia không gian

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
3

Việc tham gia hai bộ dữ liệu phi không gian dựa trên một biến ‘khóa được chia sẻ, như được mô tả trong Phần & NBSP; 3.3.3. Tham gia dữ liệu không gian áp dụng cùng một khái niệm, nhưng thay vào đó dựa vào các mối quan hệ không gian, được mô tả trong phần trước. Như với dữ liệu thuộc tính, tham gia thêm các cột mới vào đối tượng đích (đối số x trong các hàm tham gia), từ một đối tượng nguồn (y).
0Quá trình này được minh họa bằng ví dụ sau: Hãy tưởng tượng bạn có mười điểm được phân phối ngẫu nhiên trên bề mặt Trái đất và bạn hỏi, cho những điểm trên đất liền, họ đang ở trong nước nào? Thực hiện ý tưởng này trong một ví dụ có thể tái tạo sẽ xây dựng các kỹ năng xử lý dữ liệu địa lý của bạn và cho thấy cách thức hoạt động của không gian. Điểm khởi đầu là tạo ra các điểm nằm ngẫu nhiên trên bề mặt Trái đất:
1hình học
2Điểm (-115.10291 36.78178)
3Điểm (-172.98891 -71.02938)
4Điểm (-13.24134 65.23272)
5Điểm (80.97621 58.85495)
6Điểm (-28.72671 -61.25002)
7Điểm (-5.24625 19.83849)
8Điểm (-175.39891 -86.34517)
9Điểm (-4.54623 -69.64082)

Điểm (159.05039 -34.99599)4.2 shows that the

0    False
1    False
2     True
dtype: bool
1 object (top left) lacks attribute data, while the world (top right) has attributes, including country names shown for a sample of countries in the legend. Spatial joins are implemented with
0    False
1    False
2     True
dtype: bool
2, as illustrated in the code chunk below. The output is the
0    False
1    False
2     True
dtype: bool
3 object which is illustrated in Figure 4.2 (bottom left). Before creating the joined dataset, we use spatial subsetting to create world_random, which contains only countries that contain random points, to verify the number of country names returned in the joined dataset should be four (see the top right panel of Figure 4.2).

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
4

iso_a2name_longĐiểm (126.28622 -62.49509)Kịch bản được minh họa trong Hình & NBSP; 4.2 cho thấy đối tượng
0    False
1    False
2     True
dtype: bool
1 (trên cùng bên trái) thiếu dữ liệu thuộc tính, trong khi thế giới (trên cùng bên phải) có các thuộc tính, bao gồm các tên quốc gia được hiển thị cho một mẫu quốc gia trong truyền thuyết. Các kết nối không gian được thực hiện với
0    False
1    False
2     True
dtype: bool
2, như được minh họa trong đoạn mã bên dưới. Đầu ra là đối tượng
0    False
1    False
2     True
dtype: bool
3 được minh họa trong Hình & NBSP; 4.2 (phía dưới bên trái). Trước khi tạo bộ dữ liệu được tham gia, chúng tôi sử dụng tập hợp không gian để tạo World_random, chỉ chứa các quốc gia có chứa các điểm ngẫu nhiên, để xác minh số lượng tên quốc gia được trả về trong bộ dữ liệu được tham gia nên là bốn (xem bảng trên bên phải của Hình & NBSP; 4.2).
lục địa...Việc tham gia hai bộ dữ liệu phi không gian dựa trên một biến ‘khóa được chia sẻ, như được mô tả trong Phần & NBSP; 3.3.3. Tham gia dữ liệu không gian áp dụng cùng một khái niệm, nhưng thay vào đó dựa vào các mối quan hệ không gian, được mô tả trong phần trước. Như với dữ liệu thuộc tính, tham gia thêm các cột mới vào đối tượng đích (đối số x trong các hàm tham gia), từ một đối tượng nguồn (y).
4Quá trình này được minh họa bằng ví dụ sau: Hãy tưởng tượng bạn có mười điểm được phân phối ngẫu nhiên trên bề mặt Trái đất và bạn hỏi, cho những điểm trên đất liền, họ đang ở trong nước nào? Thực hiện ý tưởng này trong một ví dụ có thể tái tạo sẽ xây dựng các kỹ năng xử lý dữ liệu địa lý của bạn và cho thấy cách thức hoạt động của không gian. Điểm khởi đầu là tạo ra các điểm nằm ngẫu nhiên trên bề mặt Trái đất:hình họcĐiểm (-115.10291 36.78178)Điểm (-172.98891 -71.02938)78.841463 51921.984639 Điểm (-13.24134 65.23272)
18Điểm (80.97621 58.85495)Điểm (-28.72671 -61.25002)Điểm (-5.24625 19.83849)Điểm (-172.98891 -71.02938)70.743659 25284.586202 Điểm (-13.24134 65.23272)
52Điểm (80.97621 58.85495)Điểm (-28.72671 -61.25002)Điểm (-5.24625 19.83849)Điểm (-172.98891 -71.02938)57.007000 1865.160622 Điểm (-13.24134 65.23272)
159Điểm (80.97621 58.85495)Điểm (-28.72671 -61.25002)Điểm (-28.72671 -61.25002)Điểm (-172.98891 -71.02938)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (80.97621 58.85495)

Điểm (-28.72671 -61.25002)

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
5

Việc tham gia hai bộ dữ liệu phi không gian dựa trên một biến ‘khóa được chia sẻ, như được mô tả trong Phần & NBSP; 3.3.3. Tham gia dữ liệu không gian áp dụng cùng một khái niệm, nhưng thay vào đó dựa vào các mối quan hệ không gian, được mô tả trong phần trước. Như với dữ liệu thuộc tính, tham gia thêm các cột mới vào đối tượng đích (đối số x trong các hàm tham gia), từ một đối tượng nguồn (y).index_rightiso_a2Kịch bản được minh họa trong Hình & NBSP; 4.2 cho thấy đối tượng
0    False
1    False
2     True
dtype: bool
1 (trên cùng bên trái) thiếu dữ liệu thuộc tính, trong khi thế giới (trên cùng bên phải) có các thuộc tính, bao gồm các tên quốc gia được hiển thị cho một mẫu quốc gia trong truyền thuyết. Các kết nối không gian được thực hiện với
0    False
1    False
2     True
dtype: bool
2, như được minh họa trong đoạn mã bên dưới. Đầu ra là đối tượng
0    False
1    False
2     True
dtype: bool
3 được minh họa trong Hình & NBSP; 4.2 (phía dưới bên trái). Trước khi tạo bộ dữ liệu được tham gia, chúng tôi sử dụng tập hợp không gian để tạo World_random, chỉ chứa các quốc gia có chứa các điểm ngẫu nhiên, để xác minh số lượng tên quốc gia được trả về trong bộ dữ liệu được tham gia nên là bốn (xem bảng trên bên phải của Hình & NBSP; 4.2).
lục địalục địa...
0Quá trình này được minh họa bằng ví dụ sau: Hãy tưởng tượng bạn có mười điểm được phân phối ngẫu nhiên trên bề mặt Trái đất và bạn hỏi, cho những điểm trên đất liền, họ đang ở trong nước nào? Thực hiện ý tưởng này trong một ví dụ có thể tái tạo sẽ xây dựng các kỹ năng xử lý dữ liệu địa lý của bạn và cho thấy cách thức hoạt động của không gian. Điểm khởi đầu là tạo ra các điểm nằm ngẫu nhiên trên bề mặt Trái đất:4.0 Quá trình này được minh họa bằng ví dụ sau: Hãy tưởng tượng bạn có mười điểm được phân phối ngẫu nhiên trên bề mặt Trái đất và bạn hỏi, cho những điểm trên đất liền, họ đang ở trong nước nào? Thực hiện ý tưởng này trong một ví dụ có thể tái tạo sẽ xây dựng các kỹ năng xử lý dữ liệu địa lý của bạn và cho thấy cách thức hoạt động của không gian. Điểm khởi đầu là tạo ra các điểm nằm ngẫu nhiên trên bề mặt Trái đất:Điểm (-172.98891 -71.02938)318622525.0 78.841463 51921.984639
1hình họcĐiểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-172.98891 -71.02938)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)
2Điểm (-115.10291 36.78178)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-172.98891 -71.02938)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)
3Điểm (-172.98891 -71.02938)18.0 Điểm (80.97621 58.85495)Điểm (-172.98891 -71.02938)143819666.0 70.743659 25284.586202
4Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-172.98891 -71.02938)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)
5Điểm (80.97621 58.85495)52.0 Điểm (80.97621 58.85495)Điểm (-172.98891 -71.02938)16962846.0 57.007000 1865.160622
6Điểm (-28.72671 -61.25002)159.0 Điểm (80.97621 58.85495)Điểm (-172.98891 -71.02938)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)
7Điểm (-5.24625 19.83849)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-172.98891 -71.02938)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)
8Điểm (-175.39891 -86.34517)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-172.98891 -71.02938)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)
9Điểm (-4.54623 -69.64082)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-172.98891 -71.02938)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)Điểm (-13.24134 65.23272)

Điểm (80.97621 58.85495)

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
6

Hướng dẫn geocomputation with python - tính toán địa lý với python

Điểm (-28.72671 -61.25002)

Điểm (-5.24625 19.83849)

Điểm (-175.39891 -86.34517)

Điểm (-4.54623 -69.64082)

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
7

t50_fidđộ caohình họcindex_rightTên
02353944 2723 Điểm (1204142.603 5049971.287)12 Southland
12354404 2820 Điểm (1234725.325 5048309.302)11 Otago
22354405 2830 Điểm (1235914.511 5048745.117)11 Otago
32369113 3033 Điểm (1235914.511 5048745.117)9 Điểm (1259701.635 5076570.049)
42362630 2749 Bờ Tây10 Điểm (1378169.600 5158491.453)
Canterbury...............
962408394 2797 ...7 Điểm (1821014.190 5647970.666)
972408395 2757 Manawatu-Wanganui7 Điểm (1821014.190 5647970.666)
982408397 2751 Manawatu-Wanganui2 Điểm (1820642.909 5648331.194)
992408406 2720 Điểm (1820659.873 5649488.235)2 Điểm (1820642.909 5648331.194)
1002408411 2732 Điểm (1820659.873 5649488.235)2 Điểm (1820642.909 5648331.194)

Điểm (1820659.873 5649488.235)

Waikato

Điểm (1822262.592 5650428.656)

Điểm (1822492.184 5650492.304)

101 hàng × 5 cột

Thứ hai, …

Tham gia các lớp không phù hợp

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
8

t50_fidđộ caohình học
642372236 3724
632372235 3717 Quan hệ khoảng cách
672372252 3688 Trong khi quan hệ tôpô là nhị phân - một tính năng giao nhau với người khác hoặc không - quan hệ khoảng cách là liên tục. Khoảng cách giữa hai đối tượng được tính toán bằng phương pháp
0    False
1    False
2     True
dtype: bool
9. Phương pháp được áp dụng trên
0     True
1    False
2    False
dtype: bool
0 (hoặc
0     True
1    False
2     True
dtype: bool
1), với đối số là một hình học
0     True
1    False
2    False
dtype: bool
2 riêng lẻ. Kết quả là
0     True
1    False
2    False
dtype: bool
3 khoảng cách cặp.

Điều này được minh họa trong đoạn mã bên dưới, tìm thấy khoảng cách giữa ba điểm cao nhất ở New Zealand:

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
9

Điểm (1369317.630 5169132.284)

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
0

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
1

Điểm (1369512.866 5168235.616)

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
2

TênĐiểm (1369381.942 5168761.875)Land_areaCanterburyMedian_incomeSex_ratiohình học
10Điểm (1378169.600 5158491.453)Canterbury44504.499091 ...30100 0.975327 ...
11OtagoCanterbury31186.309188 ...26300 0.951169 ...

Điểm (1821014.190 5647970.666)

Manawatu-Wanganui

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
3

1011
0123537.158269 15497.717252
194282.773074 0.000000
293018.560814 0.000000

Điểm (1820642.909 5648331.194)

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
4

Hướng dẫn geocomputation with python - tính toán địa lý với python

Điểm (1820659.873 5649488.235)

WaikatoSection 3.4) demonstrated how to retrieve values associated with specific cell IDs or row and column combinations. Raster objects can also be extracted by location (coordinates) and other spatial objects. To use coordinates for subsetting, we can use the

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
07 method of a
from shapely.geometry import Point, LineString, Polygon
points = gpd.GeoSeries([Point(0.2,0.1), Point(0.7,0.2), Point(0.4,0.8)])
line = LineString([(0.4,0.2), (1,0.5)])
poly = Polygon([(0,0), (0,1), (1,1), (1,0.5), (0,0)])

p = gpd.GeoSeries(poly).plot(color='yellow')
p = gpd.GeoSeries(line).plot(ax=p, color='red')
points.plot(ax=p)
5 file connection object, combined with a list of coordinate tuples. The methods is demonstrated below to find the value of the cell that covers a point located at coordinates of 0.1, 0.1 in
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
09. The returned object is a generator:

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
5

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
6

Điểm (1822262.592 5650428.656)

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
7

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
8

Điểm (1822492.184 5650492.304)

Điểm (1822492.184 5650492.304)

101 hàng × 5 cột

Thứ hai, …

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
0

Tham gia các lớp không phù hợp

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
1

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
2

Section 3.4.2):

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
3

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
4

Kết quả được hiển thị trong Hình & NBSP; 4.3.4.3.

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
5

Hướng dẫn geocomputation with python - tính toán địa lý với python

Cách tiếp cận trên cũng có thể được sử dụng để thay thế một số giá trị (ví dụ: dự kiến ​​sẽ sai) bằng

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
23:

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
6

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
7

Các hoạt động này trên thực tế là các hoạt động cục bộ Boolean, vì chúng tôi so sánh hai rasters thông minh. Tiểu mục tiếp theo khám phá những hoạt động này và các hoạt động liên quan chi tiết hơn.

Bản đồ đại số

Thuật ngữ ‘Bản đồ đại số, được đặt ra vào cuối những năm 1970 để mô tả một tập hợp các quy ước, khả năng và kỹ thuật của người Hồi giáo để phân tích dữ liệu vector địa lý và (mặc dù ít nổi bật hơn) (Tomlin 1994). Trong bối cảnh này, chúng tôi xác định đại số bản đồ hẹp hơn, là các hoạt động sửa đổi hoặc tóm tắt các giá trị tế bào raster, với tham chiếu đến các tế bào xung quanh, khu vực hoặc chức năng thống kê áp dụng cho mọi ô.

Các hoạt động đại số bản đồ có xu hướng nhanh, bởi vì các bộ dữ liệu raster chỉ lưu trữ tọa độ, do đó, câu ngạn ngữ cũ raster nhanh hơn nhưng vector là người chỉnh sửa. Vị trí của các ô trong các bộ dữ liệu raster có thể được tính toán bằng cách sử dụng vị trí ma trận của nó và độ phân giải và nguồn gốc của bộ dữ liệu (được lưu trữ trong tiêu đề). Tuy nhiên, đối với việc xử lý, vị trí địa lý của một ô hầu như không phù hợp miễn là chúng tôi đảm bảo rằng vị trí tế bào vẫn giống nhau sau khi xử lý. Ngoài ra, nếu hai hoặc nhiều bộ dữ liệu raster chia sẻ cùng một mức độ, phép chiếu và độ phân giải, người ta có thể coi chúng là ma trận để xử lý.

Đây là cách mà bản đồ đại số hoạt động với gói Terra. Đầu tiên, các tiêu đề của các bộ dữ liệu raster được truy vấn và (trong trường hợp các hoạt động đại số MAP hoạt động trên nhiều bộ dữ liệu) được kiểm tra để đảm bảo các bộ dữ liệu tương thích. Thứ hai, đại số bản đồ giữ lại cái gọi là thư từ định vị một-một, có nghĩa là các tế bào không thể di chuyển. Điều này khác với đại số ma trận, trong đó các giá trị thay đổi vị trí, ví dụ khi nhân hoặc chia ma trận.

Bản đồ đại số (hoặc mô hình bản đồ với dữ liệu raster) chia các hoạt động raster thành bốn lớp con (Tomlin 1990), với mỗi lần làm việc trên một hoặc một số lưới đồng thời:

  • Hoạt động cục bộ hoặc trên mỗi tế bào
  • Các hoạt động đầu mối hoặc khu phố. Thông thường, giá trị ô đầu ra là kết quả của khối ô đầu vào 3 x 3
  • Các hoạt động Zonal tương tự như các hoạt động tiêu cự, nhưng lưới pixel xung quanh mà các giá trị mới được tính toán có thể có kích thước và hình dạng không đều
  • Hoạt động toàn cầu hoặc trên mỗi rệp; Điều đó có nghĩa là tế bào đầu ra có được giá trị của nó có khả năng từ một hoặc một số rasters

Kiểu chữ này phân loại các hoạt động đại số bản đồ theo số lượng ô được sử dụng cho từng bước xử lý pixel và loại đầu ra. Để hoàn thiện, chúng ta nên đề cập rằng các hoạt động raster cũng có thể được phân loại theo kỷ luật như địa hình, phân tích thủy văn hoặc phân loại hình ảnh. Các phần sau đây giải thích làm thế nào mỗi loại hoạt động đại số bản đồ có thể được sử dụng, có liên quan đến các ví dụ làm việc.

Hoạt động địa phương

Các hoạt động cục bộ bao gồm tất cả các hoạt động của tế bào trong một hoặc một số lớp. Đại số Raster là trường hợp sử dụng cổ điển của các hoạt động cục bộ - điều này bao gồm thêm hoặc trừ các giá trị từ một raster raster, bình phương và đa nhân. Đại số raster cũng cho phép các hoạt động logic như tìm tất cả các ô raster lớn hơn một giá trị cụ thể (5 trong ví dụ của chúng tôi dưới đây). Các hoạt động cục bộ được áp dụng bằng cú pháp hoạt động mảng

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
24, như được trình bày dưới đây:

Đầu tiên, chúng ta cần đọc các giá trị raster:

canterbury = nz[nz["Name"] == "Canterbury"]
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
canterbury_height = nz_height[sel]
9

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
0

Bây giờ, bất kỳ hoạt động mảng nguyên tố khôn ngoan có thể được áp dụng. Ví dụ:

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
0

Hình & NBSP; 4.4 thể hiện thêm một vài ví dụ.4.4 demonstrates a few more examples.

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
1

Hướng dẫn geocomputation with python - tính toán địa lý với python

Một ví dụ điển hình khác về các hoạt động cục bộ là phân loại các khoảng của các giá trị số thành các nhóm như nhóm mô hình độ cao kỹ thuật số thành thấp (lớp 1), giữa (lớp 2) và độ cao cao (lớp 3). Ở đây, chúng tôi gán các giá trị raster trong phạm vi 0 Ném12, 12 Ném24 và 24 Tiết36 được phân loại lại để lấy các giá trị 1, 2 và 3, tương ứng.

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
2

Kết quả được phân loại lại được hiển thị trong Hình & NBSP; 4.5.4.5.

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
3

Hướng dẫn geocomputation with python - tính toán địa lý với python

Việc tính toán chỉ số thực vật khác biệt được chuẩn hóa (NDVI) là hoạt động raster raster địa phương nổi tiếng (pixel-by-pixel). Nó trả về một raster với các giá trị từ -1 đến 1; Giá trị dương cho thấy sự hiện diện của thực vật sống (chủ yếu> 0,2). NDVI được tính toán từ các dải hình ảnh được cảm nhận từ xa và hồng ngoại (NIR) từ xa, thường là từ các hệ thống vệ tinh như Landsat hoặc Sentinel. Thảm thực vật hấp thụ ánh sáng rất nhiều trong phổ ánh sáng nhìn thấy, và đặc biệt là trong kênh màu đỏ, trong khi phản chiếu ánh sáng NIR, giải thích công thức NVDI:

\ [Ndvi = \ frac {NIR-red} {NIR+Red} \]

Hãy để tính toán NDVI cho tập tin vệ tinh đa năng của Công viên quốc gia Zion.

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
4

Chuyển đổi các giá trị> 1 thành dữ liệu không có dữ liệu:

Khi vẽ hình ảnh RGB bằng hàm

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
25, hàm giả định rằng:

  • Các giá trị nằm trong phạm vi
    nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    26 cho phao hoặc
    nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    27 cho số nguyên (khác được cắt)
  • Thứ tự của các ban nhạc là RGB

Để chuẩn bị cho raster nhiều băng tần cho

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
25, do đó chúng tôi đảo ngược thứ tự của các dải (ban đầu là BGR+NIR) và chia cho tối đa để đặt giá trị tối đa ở mức
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
29:

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
5

Kết quả được hiển thị trong Hình & NBSP; 4.6.4.6.

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
6

Hướng dẫn geocomputation with python - tính toán địa lý với python

Hoạt động đầu mối

Trong khi các chức năng cục bộ hoạt động trên một ô, mặc dù có thể từ nhiều lớp, các hoạt động tiêu cự có tính đến một tế bào trung tâm (tiêu cự) và hàng xóm của nó. Vùng lân cận (cũng được đặt tên là kernel, bộ lọc hoặc cửa sổ di chuyển) đang được xem xét thường là có kích thước 3 x 3 (đó là tế bào trung tâm và tám hàng xóm xung quanh của nó), nhưng có thể đảm nhận bất kỳ hình chữ nhật nào khác (không nhất thiết được xác định bởi người dùng. Hoạt động tiêu cự áp dụng hàm tổng hợp cho tất cả các ô trong vùng lân cận được chỉ định, sử dụng đầu ra tương ứng làm giá trị mới cho ô trung tâm và chuyển sang ô trung tâm tiếp theo (Hình). Các tên khác cho hoạt động này là lọc không gian và tích chập (Burrough, McDonnell và Lloyd 2015).

Trong Python, gói

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
30 có một bộ sưu tập toàn diện các chức năng để thực hiện lọc các mảng
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
24, chẳng hạn như:

  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    32
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    33
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    34 (tức là bộ lọc trung bình)
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    35 v.v ...

Trong nhóm các chức năng này, chúng tôi xác định hình dạng của cửa sổ di chuyển với một trong hai:

  • ________ 136 - một số hoặc tuple duy nhất, ngụ ý một bộ lọc có kích thước đó
  • ________ 137 - một mảng boolean, đại diện cho cả hình dạng cửa sổ và danh tính của các yếu tố được đưa vào

Ngoài các bộ lọc tích hợp cụ thể,

  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    38 Áp dụng hàm tổng sau khi nhân với mảng
    nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    39 tùy chỉnh
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    40 giúp có thể vượt qua bất kỳ chức năng tùy chỉnh nào, trong đó người dùng có thể chỉ định bất kỳ loại tính toán dựa trên cửa sổ tùy chỉnh nào.

Ví dụ: ở đây chúng tôi áp dụng bộ lọc tối thiểu với kích thước cửa sổ là

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
41 trên
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
09:

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
0

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
8

non_canterbury_height = nz_height.overlay(canterbury, how='difference')
9

Chăm sóc đặc biệt nên được trao cho các pixel cạnh. Làm thế nào họ nên được tính toán?

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
30 cung cấp một số tùy chọn thông qua tham số
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
44:

  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    45 (mặc định)
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    46
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    47
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    48
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    49

Đôi khi các cạnh raster mở rộng một cách giả tạo được coi là không phù hợp. Nói cách khác, chúng ta có thể muốn raster kết quả chứa các giá trị pixel chỉ với các cửa sổ hoàn chỉnh của Google, ví dụ như có kích thước mẫu đồng nhất hoặc vì các giá trị theo mọi hướng có vấn đề (chẳng hạn như trong tính toán địa hình). Không có tùy chọn cụ thể để không mở rộng các cạnh trong

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
30. Tuy nhiên, để có được hiệu ứng tương tự, các cạnh của mảng được lọc có thể được gán bằng
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
23, trong một số hàng và cột theo kích thước bộ lọc. Ví dụ: khi sử dụng bộ lọc
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
52, lớp pixel đầu tiên có thể được gán bằng
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
23, phản ánh thực tế là các pixel này có 3*3 vùng lân cận:

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
0

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
1

Chúng tôi có thể nhanh chóng kiểm tra xem đầu ra có đáp ứng mong đợi của chúng tôi không. Trong ví dụ của chúng tôi, giá trị tối thiểu phải luôn là góc trên bên trái của cửa sổ di chuyển (hãy nhớ rằng chúng tôi đã tạo ra raster đầu vào bằng cách tăng các giá trị ô bằng một bắt đầu ở góc trên bên trái).

Các chức năng tiêu cự hoặc bộ lọc đóng vai trò chi phối trong xử lý hình ảnh. Bộ lọc thông thấp hoặc làm mịn sử dụng hàm trung bình để loại bỏ các thái cực. Trong trường hợp dữ liệu phân loại, chúng ta có thể thay thế giá trị trung bình bằng chế độ, là giá trị phổ biến nhất. Ngược lại, các bộ lọc thông cao làm nổi bật các tính năng. Các bộ lọc Laplace và Sobel phát hiện dòng có thể đóng vai trò là một ví dụ ở đây.

Xử lý địa hình, tính toán các đặc điểm địa hình như độ dốc, hướng và hướng dòng chảy, dựa vào các chức năng tiêu cự. Hàm

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
54 từ gói
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
55 có thể được sử dụng để tính toán các số liệu chung, được chỉ định thông qua đối số
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
56, cụ thể là:

  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    57 Horn (1981) DOI: 10.1109/Proc.1981.11918
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    58 Horn (1981) DOI: 10.1109/Proc.1981.11918
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    59 Horn (1981) DOI: 10.1109/Proc.1981.11918
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    60 Horn (1981) DOI: 10.1109/Proc.1981.11918
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    61 Horn (1981) DOI: 10.1109/Proc.1981.11918
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    62 Zevenbergen và Thorne (1987) DOI: 10.1002/ESP.3290120107
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    63 Zevenbergen và Thorne (1987) DOI: 10.1002/ESP.3290120107
  • nz = gpd.read_file("data/nz.gpkg")
    nz_height = gpd.read_file("data/nz_height.gpkg")
    world = gpd.read_file('data/world.gpkg')
    src_elev = rasterio.open("data/elev.tif")
    src_multi_rast = rasterio.open("data/landsat.tif")
    src_grain = rasterio.open('data/grain.tif')
    64 Zevenbergen và Thorne (1987) DOI: 10.1002/ESP.3290120107

Hoạt động khu vực

Giống như các hoạt động đầu mối, các hoạt động khu vực áp dụng chức năng tổng hợp cho nhiều ô raster. Tuy nhiên, raster thứ hai, thường là với các giá trị phân loại, xác định các bộ lọc khu vực (hoặc ‘vùng,) trong trường hợp hoạt động khu vực, trái ngược với cửa sổ lân cận được xác định trước trong trường hợp hoạt động tiêu cự được trình bày trong phần trước. Do đó, các tế bào raster xác định bộ lọc khu vực không nhất thiết phải là hàng xóm. Raster kích thước hạt của chúng tôi là một ví dụ điển hình, như được minh họa trong bảng bên phải của Hình 3.2: Kích thước hạt khác nhau được lan truyền không đều trong suốt raster. Cuối cùng, kết quả của một hoạt động khu vực là một bảng tóm tắt được nhóm theo khu vực, đó là lý do tại sao hoạt động này còn được gọi là thống kê khu vực trong thế giới GIS. Điều này trái ngược với các hoạt động tiêu cự trả về một đối tượng raster.

Để chứng minh, chúng ta hãy quay lại các raster

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
65 và
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
09 (Hình 3.2). Để tính toán số liệu thống kê khu vực, chúng tôi sử dụng các mảng có giá trị raster. Mảng
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
09 đã được nhập sớm hơn:

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[0])
nz_height.plot(ax=base, color="None", edgecolor="red")
base = nz.plot(color="white", edgecolor="lightgrey", ax=axes[1])
canterbury_height.plot(ax=base, color="None", edgecolor="red");
0

Mảng

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
65 được nhập ở đây bằng cách sử dụng
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
69:

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
3

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
4

Việc quốc tế của chúng tôi là tính toán mức trung bình (hoặc bất kỳ hàm tóm tắt nào khác, cho vấn đề đó) của độ cao trong mỗi vùng được xác định bởi các giá trị hạt. Đầu tiên, chúng ta có thể thu được các giá trị duy nhất xác định các khu vực bằng cách sử dụng

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
70:

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
5

Bây giờ, chúng ta có thể sử dụng sự hiểu biết từ điển để chia ra mảng

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
09 thành các mảng một chiều riêng biệt với các giá trị trên một nhóm
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
65, với các khóa là các giá trị
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
65 duy nhất:

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
6

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
7

Ở giai đoạn này, chúng ta có thể mở rộng biểu thức hiểu từ điển để tính độ cao trung bình liên quan đến từng lớp kích thước hạt. Thay vì đặt các giá trị độ cao (

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
74) vào các giá trị từ điển, chúng tôi đặt giá trị trung bình của chúng (
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
75):

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
8

base = nz.plot(color="white", edgecolor="lightgrey")
non_canterbury_height.plot(ax=base, color="None", edgecolor="red");
9

Điều này trả về các số liệu thống kê cho từng loại, ở đây độ cao trung bình cho mỗi lớp kích thước hạt. Ví dụ, độ cao trung bình trong các pixel được đặc trưng bởi kích thước hạt

nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
76 là
nz = gpd.read_file("data/nz.gpkg")
nz_height = gpd.read_file("data/nz_height.gpkg")
world = gpd.read_file('data/world.gpkg')
src_elev = rasterio.open("data/elev.tif")
src_multi_rast = rasterio.open("data/landsat.tif")
src_grain = rasterio.open('data/grain.tif')
77, v.v.

Hoạt động và khoảng cách toàn cầu

Bản đồ các đối tác đại số trong xử lý vector

Bản đồ các đối tác đại số trong xử lý vector

Bản đồ các đối tác đại số trong xử lý vector

  • Hợp nhất Rasters