Geocomputation （3）Spatial data operations

news2024/11/13 15:04:23

Geocomputation （3）Spatial data operations

来源：https://github.com/geocompx/geocompy

1.准备

#| echo: false
import pandas as pd
import matplotlib.pyplot as plt
pd.set_option("display.max_rows", 4)
pd.set_option("display.max_columns", 6)
pd.options.display.max_colwidth = 35
plt.rcParams["figure.figsize"] = (5, 5)

import shapely.geometry
import geopandas as gpd
import numpy as np
import os
import rasterio
import scipy.ndimage
import rasterio.merge

from rasterio.plot import show

读取数据：

#| echo: false
from pathlib import Path
data_path = Path("data")
file_path = Path("data/landsat.tif")
if not file_path.exists():
  if not data_path.is_dir():
     os.mkdir(data_path)
  print("Attempting to get the data")
  import requests
  r = requests.get("https://github.com/geocompx/geocompy/releases/download/0.1/landsat.tif")  
  with open(file_path, "wb") as f:
    f.write(r.content)

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

2.简介

3.矢量数据的空间操作

3.1 空间子集化

空间子集化是获取空间对象并返回仅包含与另一个对象在空间中相关的特征的新对象的过程。
类似于属性子集（在 @sec-vector-attribute-subsetting 中介绍），可以使用方括号（“[”）运算符使用语法“x[y]”创建“GeoDataFrame”的子集，其中“x”是 一个“GeoDataFrame”，将从中返回行/特征的子集，“y”是“子集对象”。
反过来，可以使用一种二元几何关系方法创建“y”，例如“.intersects”（请参阅@sec-topological-relations）。

为了演示空间子集化，我们将使用投影坐标系中的“nz”和“nz_height”图层，它们分别包含新西兰 16 个主要地区和 101 个最高点的地理数据 (@fig-spatial-subset) .
以下代码行创建一个表示坎特伯雷 (canterbury) 的对象，然后使用空间子集返回该区域中的所有高点 (canterbury_height)：

canterbury = nz[nz['Name'] == 'Canterbury']
canterbury
NameIslandLand_area...Median_incomeSex_ratiogeometry
10CanterburySouth44504.499091...301000.975327MULTIPOLYGON (((1686901.914 535...

1 rows × 7 columns

# Does each 'nz_height' point intersect with 'canterbury'?
sel = nz_height.intersects(canterbury['geometry'].iloc[0])
sel

0      False
1      False
       ...  
99     False
100    False
Length: 101, dtype: bool

canterbury_height = nz_height[sel]
canterbury_height
t50_fidelevationgeometry
423626302749POINT (1378169.600 5158491.453)
523628142822POINT (1389460.041 5168749.086)
............
9323803002711POINT (1654213.379 5349962.973)
9423803082885POINT (1654898.622 5350462.779)

70 rows × 3 columns

#| label: fig-spatial-subset
#| fig-cap: Spatial subsetting of points by intersection with polygon

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.plot(ax=base, color='grey', edgecolor='none')
canterbury_height.plot(ax=base, color='None', edgecolor='red')
axes[0].set_title('Original')
axes[1].set_title('Subset (intersects)');

在这里插入图片描述

与属性子集 (@sec-vector-attribute-subsetting) 一样，我们使用布尔系列 (sel)，其长度与过滤表 (nz_height) 中的行数相同，基于 条件应用于自身。
不同的是，条件不是属性值的比较，而是空间关系的评估。
也就是说，我们使用“.intersects”方法评估“nz_height”的每个几何图形是否与“canterbury”几何图形相交。

各种拓扑关系可用于空间子集化，确定目标对象中的要素与要选择的子集化对象必须具有的空间关系类型。
这些包括接触、交叉或内部，我们很快就会在@sec-topological-relations 中看到。
我们在上一个示例中使用的最常用的方法 .intersects 方法是“捕获所有”拓扑关系，它将返回目标中接触、交叉或位于源“子集”对象内的特征。
或者，我们可以评估其他方法，例如“.disjoint”来获取所有“不”与坎特伯雷相交的点：

# Is each 'nz_height' point disjoint from 'canterbury'?
sel = nz_height.disjoint(canterbury['geometry'].iloc[0])
canterbury_height2 = nz_height[sel]

如@fig-spatial-subset-disjoint 所示：

#| label: fig-spatial-subset-disjoint
#| fig-cap: Spatial subsetting of points disjoint from a polygon

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.plot(ax=base, color='grey', edgecolor='none')
canterbury_height2.plot(ax=base, color='None', edgecolor='red');
axes[0].set_title('Original')
axes[1].set_title('Subset (disjoint)');

在这里插入图片描述

如果我们需要同时根据多个几何图形进行子集化，例如找出哪些点与坎特伯雷和南地相交，我们可以在应用“.intersects”（或任何其他）运算符之前溶解过滤子集：

canterbury_southland = nz[(nz['Name'] == 'Canterbury') | (nz['Name'] == 'Southland')]
canterbury_southland = canterbury_southland.unary_union
canterbury_southland

在这里插入图片描述

sel = nz_height.intersects(canterbury_southland)
canterbury_southland_height = nz_height[sel]
canterbury_southland_height
t50_fidelevationgeometry
023539442723POINT (1204142.603 5049971.287)
423626302749POINT (1378169.600 5158491.453)
............
9323803002711POINT (1654213.379 5349962.973)
9423803082885POINT (1654898.622 5350462.779)

71 rows × 3 columns

生成的子集显示在 @fig-spatial-subset2 中：

#| label: fig-spatial-subset2
#| fig-cap: Spatial subsetting of points by intersection with more that one polygon

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])
gpd.GeoSeries(canterbury_southland).plot(ax=base, color='grey', edgecolor='none')
canterbury_southland_height.plot(ax=base, color='None', edgecolor='red')
axes[0].set_title('Original')
axes[1].set_title('Subset (intersects)');

在这里插入图片描述

下一节将进一步探讨不同类型的空间关系，也称为二元谓词（其中“.intersects”和“.disjoint”是两个示例），可用于识别两个要素是否在空间上相关。

3.2 拓扑关系

拓扑关系描述了对象之间的空间关系。 “二元拓扑关系”，全称为“二元拓扑关系”，是关于由有序点集（通常形成点、 线和多边形）在两个或多个维度（Egenhofer 和 Herring 1990）。 这听起来可能相当抽象，事实上，拓扑关系的定义和分类是基于 1966 年首次以书籍形式出版的数学基础（Spanier 1995），代数拓扑领域一直延续到 21 世纪（Dieck 2008）。

尽管拓扑关系起源于数学，但可以参考测试常见类型空间关系的常用函数的可视化来直观地理解拓扑关系。 图 4.2 显示了各种几何对及其关联关系。 图 4.2 中的第三对和第四对（从左到右，然后向下）表明，对于某些关系，顺序很重要：而相等、相交、交叉、接触和重叠的关系是对称的，这意味着如果 function(x, y) 为 true，则 function(y, x) 也为 true，其中几何顺序很重要的关系（例如 contains 和 inside 不重要）。 请注意，每个几何对都有一个“DE-9IM”字符串，例如 FF2F11212，如下一节所述。

矢量几何之间的拓扑关系，受到 Egenhofer 和 Herring (1990) 中图 1 和图 2 的启发。 为每个几何对打印函数 (x, y) 为真的关系，其中 x 用粉色表示，y 用蓝色表示。 每对空间关系的性质由维度扩展 9 相交模型字符串描述。

在shapely中，针对不同类型的拓扑关系进行测试的函数被称为“关系”。 为了了解拓扑关系在实践中如何工作，让我们创建一个简单的可重现示例，以图 4.2 中所示的关系为基础，并巩固上一章中关于矢量几何如何表示的知识（@sec-geometry-columns 和 @sec-geometries） :

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

生成的名为“points”、“line”和“poly”的“GeoSeries”可视化如下：

#| label: fig-spatial-relations-geoms
#| fig-cap: Points, line and polygon objects arranged to illustrate topological relations

ax = gpd.GeoSeries(poly).plot(color='grey', edgecolor='red')
ax = gpd.GeoSeries(line).plot(ax=ax, color='black')
points.plot(ax=ax, color='none', edgecolor='black')
for x, y, label in zip(points.x, points.y, [1,2,3]):
    ax.annotate(label, xy=(x, y), xytext=(3, 3), weight='bold', textcoords="offset points")

在这里插入图片描述

一个简单的查询是：“points”中的哪些点以某种方式与多边形“poly”相交？ 这个问题可以通过检查来回答（点 1 和点 3 分别接触并在多边形内）。 这个问题可以用.intersects方法来回答，如下。 请注意，我们评估“GeoSeries”（“点”）中的每个几何图形与单个“shapely”几何图形（“poly.iloc[0]”）之间的关系。 当两个输入都是“GeoSeries”时，在按索引（“align=True”，默认值）或按位置（“align=False”）对齐的几何图形之间进行成对评估：

points.intersects(poly.iloc[0])

0     True
1    False
2     True
dtype: bool

结果是一个布尔值“Series”。 它的内容应该符合您的直觉：第一个和第三个点返回正（“True”）结果，第二个点返回负结果（“False”）在多边形边界之外。

同样，上面的输出是一个“Series”，其中每个值代表第一个输入（“points”）中的一个特征。 如果我们需要获得成对结果的矩阵，我们可以使用“.apply”，如下所示。 在这种情况下，结果是一个“DataFrame”，其中每行代表一个“points”几何图形，每列代表一个“poly”几何图形（在本例中只有一个）L

points.apply(lambda x: poly.intersects(x))
0
0True
1False
2True

即使在特征刚好接触的情况下，“.intersects”也会返回“True”：intersects 是一种“包罗万象”的拓扑运算，可识别多种类型的空间关系，如图 4.2 所示。 更具限制性的问题包括哪些点位于多边形内，哪些要素位于多边形上或包含与其共享的边界？ 这些可以回答如下：

points.within(poly.iloc[0])

0    False
1    False
2     True
dtype: bool

points.touches(poly.iloc[0])

0     True
1    False
2    False
dtype: bool

请注意，虽然第一个点接触到边界多边形，但它不在其中； 第三个点在多边形内但不接触其边界的任何部分。 .intersects 的对立面是 .disjoint，它只返回在空间上与选择对象没有任何关系的对象：

points.disjoint(poly.iloc[0])

0    False
1     True
2    False
dtype: bool

另一种有用的关系类型是“在距离内”，我们在其中检测与缓冲特定距离的目标相交的特征。 缓冲距离决定了目标对象在被选中之前需要多近。 这可以通过字面上缓冲（@sec-geometries）目标几何，并评估交集（.intersects，见上文）来完成。 另一种方法是计算距离并将它们与距离阈值进行比较：

points.distance(poly.iloc[0]) < 0.2

0    True
1    True
2    True
dtype: bool

请注意，虽然点 2 距离“poly”最近的顶点超过 0.2 个单位的距离，但当距离设置为 0.2 时，它仍然被选中。 这是因为距离是测量到最近的边缘，在本例中是位于图 4.2 中点 2 正上方的多边形部分。 我们可以验证点2到多边形的实际距离是0.13，如下：

points.iloc[1].distance(poly.iloc[0])

0.13416407864998736

3.3 空间连接

连接两个非空间数据集依赖于共享的“键”变量，如 @sec-vector-attribute-joining 中所述。 空间数据连接应用相同的概念，但依赖于空间关系，如上一节所述。 与属性数据一样，连接将新列从源对象 (y) 添加到目标对象（连接函数中的参数 x）。

下面的例子说明了这个过程：假设你有十个随机分布在地球表面的点，你问，对于陆地上的点，它们位于哪些国家？ 在可重复的示例中实现这个想法将培养您的地理数据处理技能并展示空间连接的工作原理。 起点是创建随机散布在地球表面上的点：

np.random.seed(11)  ## set seed for reproducibility
bb = world.total_bounds  ## the world's bounds
x = np.random.uniform(low=bb[0], high=bb[2], size=10)
y = np.random.uniform(low=bb[1], high=bb[3], size=10)
random_points = gpd.points_from_xy(x, y, crs=4326)
random_points = gpd.GeoSeries(random_points)
random_points = gpd.GeoDataFrame({'geometry': random_points})
random_points
geometry
0POINT (-115.10291 36.78178)
1POINT (-172.98891 -71.02938)
......
8POINT (159.05039 -34.99599)
9POINT (126.28622 -62.49509)

10 rows × 1 columns

@fig-spatial-join 中所示的场景显示“random_points”对象（左上）缺少属性数据，而世界（右上）具有属性，包括图例中国家/地区样本显示的国家/地区名称。 空间连接是通过“gpd.sjoin”实现的，如下面的代码块所示。 输出是“random_joined”对象，如@fig-spatial-join（左下）所示。 在创建连接数据集之前，我们使用空间子集创建 world_random，其中仅包含包含随机点的国家/地区，以验证连接数据集中返回的国家/地区名称数量应为 4（请参阅 @fig-spatial- 的右上角面板） 加入）。

# Subset
world_random = world[world.intersects(random_points.unary_union)]
world_random
iso_a2name_longcontinent...lifeExpgdpPercapgeometry
4USUnited StatesNorth America...78.84146351921.984639MULTIPOLYGON (((-171.73166 63.7...
18RURussian FederationEurope...70.74365925284.586202MULTIPOLYGON (((-180.00000 64.9...
52MLMaliAfrica...57.0070001865.160622MULTIPOLYGON (((-11.51394 12.44...
159AQAntarcticaAntarctica...NaNNaNMULTIPOLYGON (((-180.00000 -89....

4 rows × 11 columns

# Spatial join
random_joined = gpd.sjoin(random_points, world, how='left')
random_joined
geometryindex_rightiso_a2...poplifeExpgdpPercap
0POINT (-115.10291 36.78178)4.0US...318622525.078.84146351921.984639
1POINT (-172.98891 -71.02938)NaNNaN...NaNNaNNaN
........................
8POINT (159.05039 -34.99599)NaNNaN...NaNNaNNaN
9POINT (126.28622 -62.49509)NaNNaN...NaNNaNNaN

10 rows × 12 columns

#| label: fig-spatial-join
#| fig-cap: Illustration of a spatial join. A new attribute variable is added to random points (top left) from source world object (top right) resulting in the data represented in the final panel.

fig, axes = plt.subplots(2, 2, figsize=(9, 5))
base = world.plot(color='white', edgecolor='lightgrey', ax=axes[0][0])
random_points.plot(ax=base, color='None', edgecolor='red')
base = world.plot(color='white', edgecolor='lightgrey', ax=axes[0][1])
world_random.plot(ax=base, column='name_long')
base = world.plot(color='white', edgecolor='lightgrey', ax=axes[1][0])
random_joined.geometry.plot(ax=base, color='grey');
random_joined.plot(ax=base, column='name_long', legend=True)
fig.delaxes(axes[1][1]);

在这里插入图片描述

3.4 非重叠连接

有时，两个地理数据集并不接触，但仍然具有很强的地理关系。
数据集“cycle_hire”和“cycle_hire_osm”提供了一个很好的例子。
对它们进行绘图表明，它们通常密切相关，但并不接触，如 @fig-cycle-hire 所示，它是使用以下代码创建的：

#| label: fig-cycle-hire
#| fig-cap: The spatial distribution of cycle hire points in London based on official data (blue) and OpenStreetMap data (red). 

base = cycle_hire.plot(edgecolor='blue', color='none')
cycle_hire_osm.plot(ax=base, edgecolor='red', color='none');

在这里插入图片描述

我们可以通过创建“.intersects”关系的成对布尔矩阵来检查是否有任何点相同，然后评估其中的任何值是否为“True”。
请注意，“.to_numpy”方法应用于从“DataFrame”到“numpy”数组，其中“.any”给出全局而不是逐行摘要，这正是我们在本例中想要的：

m = cycle_hire['geometry'].apply(
  lambda x: cycle_hire_osm['geometry'].intersects(x)
)
m.to_numpy().any()

False

想象一下，我们需要将 cycle_hire_osm 中的容量变量连接到 cycle_hire 中包含的官方“目标”数据。
这是需要非重叠连接的时候。
空间连接 (gpd.sjoin) 以及缓冲几何可用于执行此操作。
下面使用 20 m 的阈值距离对此进行了演示。
请注意，我们将数据转换为投影 CRS (27700) 以使用以米为单位的实际缓冲距离。

crs = 27700
cycle_hire2 = cycle_hire.copy().to_crs(crs)
cycle_hire2['geometry'] = cycle_hire2.buffer(20)
z = gpd.sjoin(
  cycle_hire2, 
  cycle_hire_osm.to_crs(crs)
)
z[['name_left','name_right']]
name_leftname_right
0River StreetRiver Street
1Phillimore GardensKensington, Phillimore Gardens
.........
729Ranelagh GardensRanelagh Gardens, Fulham
737Tallis StreetTallis Street, Temple

458 rows × 2 columns

结果“z”表明目标对象“cycle_hire”在“cycle_hire_osm”阈值距离内有438个点。
请注意，连接结果中的行数大于目标。
这是因为 cycle_hire 中的某些自行车租赁站在 cycle_hire_osm 中有多个匹配项。
为了聚合重叠点的值并返回平均值，我们可以使用@attr中学到的聚合方法，从而得到一个与目标具有相同行数的对象。
我们还使用“.centroid”从缓冲区返回到点：

z = z[['id', 'capacity', 'geometry']] \
    .dissolve(by='id', aggfunc='mean') \
    .reset_index()
z['geometry'] = z.centroid
z
idgeometrycapacity
01POINT (531203.517 182832.066)9.0
12POINT (525208.067 179391.922)27.0
............
436765POINT (524652.998 175817.001)29.0
437773POINT (531435.032 180916.010)14.0

438 rows × 3 columns

可以通过将源“cycle_hire_osm”数据的容量图与新对象“z”（@fig-cycle-hire-z）中的结果进行比较来验证附近站点的容量：

#| label: fig-cycle-hire-z
#| fig-cap: Non-overlapping join input (left) and result (right)

fig, axes = plt.subplots(ncols=2, figsize=(8,3))
cycle_hire_osm.plot(column='capacity', legend=True, ax=axes[0])
z.plot(column='capacity', legend=True, ax=axes[1])
axes[0].set_title('Input')
axes[1].set_title('Join result');

在这里插入图片描述

3.5 空间聚合

与属性数据聚合一样，空间数据聚合会压缩数据：聚合输出的行数少于非聚合输入的行数。
统计聚合函数（例如均值、平均值或总和）汇总一个变量的多个值，并为每个分组变量返回一个值。
@sec-vector-attribute-aggregation 演示了 .groupby 方法如何结合诸如 .sum 之类的汇总函数，基于属性变量压缩数据。
本节展示如何使用空间连接与非空间聚合相结合来实现按空间对象分组。

回到新西兰的例子，假设您想要找出每个地区高点的平均高度。
源 (nz) 的几何形状定义了目标对象 (nz_height) 中的值如何分组。
这可以通过三个步骤完成：

  1. 找出每个 nz_height 点落在哪个 nz 区域—使用 gpd.sjoin
  2. 汇总每个区域的平均海拔—使用 .groupby 和 .mean
  3. 将结果连接回 nz — 使用 pd.merge

首先，我们使用空间连接“附加”每个点的区域分类。
请注意，我们使用的是所需的最小列集：几何图形（用于空间连接）、点高程（稍后计算平均值）和区域名称（在将结果连接回 时用作键 新西兰）。

nz_height2 = gpd.sjoin(
  nz_height[['elevation', 'geometry']], 
  nz[['Name', 'geometry']], 
  how='left'
)
nz_height2
elevationgeometryindex_rightName
02723POINT (1204142.603 5049971.287)12Southland
12820POINT (1234725.325 5048309.302)11Otago
...............
992720POINT (1822262.592 5650428.656)2Waikato
1002732POINT (1822492.184 5650492.304)2Waikato

101 rows × 4 columns

其次，我们计算平均海拔：

nz_height3 = nz_height2.groupby('Name')[['elevation']].mean()
nz_height3
elevation
Name
Canterbury2994.600000
Manawatu-Wanganui2777.000000
......
Waikato2734.333333
West Coast2889.454545

7 rows × 1 columns

第三步也是最后一步是将平均值与“nz”层连接起来：

nz_height4 = pd.merge(nz[['Name', 'geometry']], nz_height3, on='Name', how='left')
nz_height4
Namegeometryelevation
0NorthlandMULTIPOLYGON (((1745493.196 600...NaN
1AucklandMULTIPOLYGON (((1803822.103 590...NaN
............
14NelsonMULTIPOLYGON (((1624866.279 541...NaN
15MarlboroughMULTIPOLYGON (((1686901.914 535...2720.0

16 rows × 3 columns

我们现在已经创建了 nz_height4 层，它给出了每个多边形的平均 nz_height 高程值。 结果显示在@fig-nz-avg-nz-height 中。 请注意，“missing_kwds”部分确定了缺少符号系统属性（“elevation”）的几何样式，因为没有“nz_height”点与它们重叠。 默认是忽略它们，这通常不是我们想要的。 使用{'color':'none','edgecolor':'black'}，这些多边形显示为黑色轮廓且无填充。

#| label: fig-nz-avg-nz-height
#| fig-cap: Average height of the top 101 high points across the regions of New Zealand

nz_height4.plot(
  column='elevation', 
  legend=True,
  missing_kwds={'color': 'none', 'edgecolor': 'black'}
);

在这里插入图片描述

3.6 连接不一致的层

空间一致性是与空间聚合相关的一个重要概念。 如果两个对象具有共享边界，则聚合对象（我们将其称为“y”）与目标对象（“x”）一致。 行政边界数据通常就是这种情况，较大的单位（例如英国的中层超级输出区 (MSOA) 或许多其他欧洲国家的地区）由许多较小的单位组成。

相比之下，不一致的聚合对象不与目标共享共同边界（Qiu、Zhang 和 Zhou 2012）。 这对于 @fig-nz-and-grid 中所示的空间聚合（和其他空间操作）来说是有问题的：聚合每个子区域的质心将不会返回准确的结果。 面积插值通过将值从一组面积单位转移到另一组面积单位来克服这个问题，使用一系列算法，包括简单的面积加权方法和更复杂的方法，例如“pycnophyctic”方法（Tobler 1979）。

为了演示，我们将创建一个“合成”层，其中包含大小为矩形的常规网格 100 × 100 100\times100 100×100 k m km km，覆盖 nz 层的范围：

xmin, ymin, xmax, ymax = nz.total_bounds
res = 100000
cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax+res)), res))
rows = list(range(int(np.floor(ymin)), int(np.ceil(ymax+res)), res))
rows.reverse()
polygons = []
for x in cols:
    for y in rows:
        polygons.append( shapely.geometry.Polygon([(x,y), (x+res, y), (x+res, y-res), (x, y-res)]) )
grid = gpd.GeoDataFrame({'geometry': polygons}, crs=nz.crs)
sel = grid.intersects(shapely.geometry.box(*nz.total_bounds))
grid = grid[sel]
grid['id'] = grid.index
grid
geometryid
0POLYGON ((1090143.000 6248536.0...0
1POLYGON ((1090143.000 6148536.0...1
.........
157POLYGON ((1990143.000 4948536.0...157
158POLYGON ((1990143.000 4848536.0...158

150 rows × 2 columns

如@fig-nz-and-grid 所示。

#| label: fig-nz-and-grid
#| fig-cap: The nz layer and a regular grid of rectangles

base = grid.plot(color='none', edgecolor='grey')
nz.plot(ax=base, column='Population', edgecolor='black', legend=True);

在这里插入图片描述

现在，我们的目标是将“人口”属性“转移”到矩形网格多边形，这是不一致层之间连接的示例。 为此，我们基本上需要计算每个“网格”单元格与该单元格重合的“nz”多边形中人口的加权总和。 加权和计算中的权重是整个“nz”多边形中重合“部分”的面积之间的比率。 也就是说，我们（不可避免地）假设每个“nz”多边形中的人口在空间上均匀分布，因此部分“nz”多边形包含相应的部分人口规模。

我们首先使用“.area”方法计算每个“nz”多边形的整个面积（参见@sec-area-length）：

nz['area'] = nz.area
nz
NameIslandLand_area...Sex_ratiogeometryarea
0NorthlandNorth12500.561149...0.942453MULTIPOLYGON (((1745493.196 600...1.289058e+10
1AucklandNorth4941.572557...0.944286MULTIPOLYGON (((1803822.103 590...4.911565e+09
........................
14NelsonSouth422.195242...0.925967MULTIPOLYGON (((1624866.279 541...4.080754e+08
15MarlboroughSouth10457.745485...0.957792MULTIPOLYGON (((1686901.914 535...1.046485e+10

16 rows × 8 columns

接下来，我们使用 .overlay 方法来计算 nz 和 grid 之间的两两交集，这里命名为 nz_grid 。 我们还计算交叉点的面积，在此命名为“area_sub”。 如果“nz”多边形完全位于单个“grid”多边形内，则“area_sub”将等于“area”； 否则，它会更小：

nz_grid = nz.overlay(grid)
nz_grid = nz_grid[['id', 'area', 'Population', 'geometry']]
nz_grid['area_sub'] = nz_grid.area
nz_grid
idareaPopulationgeometryarea_sub
0641.289058e+10175500.0POLYGON ((1590143.000 6162776.6...3.231015e+08
1801.289058e+10175500.0POLYGON ((1600430.924 6148536.0...4.612641e+08
..................
108874.080754e+0851400.0POLYGON ((1643904.326 5448536.0...1.716260e+07
109871.046485e+1046200.0MULTIPOLYGON (((1676051.777 544...4.526248e+08

110 rows × 5 columns

我们将在…中详细说明 .overlay 方法。

生成的图层“nz_grid”显示在@fig-nz-and-grid2 中。

#| label: fig-nz-and-grid2
#| fig-cap: "The nz layer and a regular grid of rectangles: after an overlay operation"

base = grid.plot(color='none', edgecolor='grey')
nz_grid.plot(ax=base, column='Population', edgecolor='black', legend=True);

在这里插入图片描述

请注意，每个“交叉点”仍然保留其“nz”“起源”特征的“Population”属性，如 @fig-nz-and-grid2 中所示。 然而，每个“nz_grid”特征的实际人口规模较小（或相等），具体取决于它在原始“nz”特征中所占的地理区域比例。 为了进行“修正”，我们首先计算比率（area_prop），然后将其乘以人口。 新的（小写）属性“population”现在可以正确估计“nz_grid”中的人口规模：

nz_grid['area_prop'] = nz_grid['area_sub'] / nz_grid['area']
nz_grid['population'] = nz_grid['Population'] * nz_grid['area_prop']
nz_grid
idareaPopulation...area_subarea_proppopulation
0641.289058e+10175500.0...3.231015e+080.0250654398.897141
1801.289058e+10175500.0...4.612641e+080.0357836279.925114
........................
108874.080754e+0851400.0...1.716260e+070.0420572161.752203
109871.046485e+1046200.0...4.526248e+080.0432521998.239223

110 rows × 7 columns

剩下要做的是总结（参见@sec-vector-attribute-aggregation）形成相同网格单元的所有部分的人口：

nz_grid = nz_grid.groupby('id')['population'].sum().reset_index()
nz_grid
idpopulation
01167.533590
11215339.996965
.........
5514931284.910446
56150129.326331

57 rows × 2 columns

并将它们连接回“网格”层（请参阅@sec-vector-attribute-joining）：

grid = pd.merge(grid, nz_grid[['id', 'population']], on='id', how='left')
grid
geometryidpopulation
0POLYGON ((1090143.000 6248536.0...0NaN
1POLYGON ((1090143.000 6148536.0...1NaN
............
148POLYGON ((1990143.000 4948536.0...157NaN
149POLYGON ((1990143.000 4848536.0...158NaN

150 rows × 3 columns

@fig-nz-and-grid3 中显示了最终结果“grid”以及来自“nz”的不一致连接的“population”属性。

#| label: fig-nz-and-grid3
#| fig-cap: "The nz layer and a regular grid of rectangles: final result"

base = grid.plot(column='population', edgecolor='black', legend=True);
nz.plot(ax=base, color='none', edgecolor='grey', legend=True);

在这里插入图片描述

我们可以证明，正如预期的那样，尽管几何形状不同，“nz”和“grid”中的总人口是相同的（因为我们创建了“grid”来完全覆盖“nz”）：

nz['Population'].sum()

4787200.0

grid['population'].sum()

4787199.999999998

本节中的过程称为空间广泛（例如人口）变量的面积加权插值。 空间密集变量（例如人口密度）的面积加权插值几乎是相同的，只是我们必须计算加权“.mean”而不是“.sum”，以保留平均值而不是总和。

3.7 距离关系

虽然拓扑关系是二元的（一个特征要么与另一个特征相交，要么不相交），但距离关系是连续的。 两个物体之间的距离是用“distance”方法计算的。 该方法应用于“GeoSeries”（或“GeoDataFrame”），参数是单个“shapely”几何图形。 结果是成对距离的“系列”。

下面的代码块对此进行了说明，该代码块查找新西兰三个最高点之间的距离：

nz_heighest = nz_height.sort_values(by='elevation', ascending=False).iloc[:3, :]
nz_heighest
t50_fidelevationgeometry
6423722363724POINT (1369317.630 5169132.284)
6323722353717POINT (1369512.866 5168235.616)
6723722523688POINT (1369381.942 5168761.875)

以及第 4.2.1 节中创建的坎特伯雷地区的地理质心：

canterbury_centroid = canterbury.centroid.iloc[0]

以下是距离：

nz_heighest.distance(canterbury_centroid)

64    115539.995747
63    115390.248038
67    115493.594066
dtype: float64

要获得距离矩阵，即对象“x”和“y”中所有特征组合之间的成对距离集，我们需要使用“.apply”方法。 这在下面的命令中进行了说明，它找到了 nz_height 中的前三个特征与对象 co 代表的新西兰奥塔哥和坎特伯雷地区之间的距离：

sel = nz['Name'].str.contains('Canter|Otag')
co = nz[sel]
co
NameIslandLand_area...Sex_ratiogeometryarea
10CanterburySouth44504.499091...0.975327MULTIPOLYGON (((1686901.914 535...4.532656e+10
11OtagoSouth31186.309188...0.951169MULTIPOLYGON (((1335204.789 512...3.190356e+10

2 rows × 8 columns

距离矩阵“d”的获取方式如下（从技术上讲，这是一个“DataFrame”）。 用简单的语言来说，我们从 nz_height.iloc[:3, :] 中获取每一行的几何形状，并以该行作为参数对 co 应用 .distance 方法：

d = nz_height.iloc[:3, :].apply(lambda x: co.distance(x['geometry']), axis=1)
d
1011
0123537.15826915497.717252
194282.7730740.000000
293018.5608140.000000

请注意，“nz_height”中的第二个和第三个特征与“co”中的第二个特征之间的距离为零。 这说明了点和多边形之间的距离是指到多边形任何部分的距离：“nz_height”中的第二个和第三个点位于奥塔哥，这可以通过绘制它们来验证：

base = co.iloc[[1]].plot(color='none')
nz_height.iloc[1:3, :].plot(ax=base);

在这里插入图片描述

4. 栅格数据的空间操作

4.1 空间子集化

上一章（@sec-manipulate-raster-objects 节）演示了如何检索与特定单元 ID 或行和列组合关联的值。 栅格对象还可以通过位置（坐标）和其他空间对象来提取。 要使用坐标进行子集化，我们可以使用“rasterio”文件连接对象的“.sample”方法，并结合坐标元组列表。 下面演示了这些方法，用于查找覆盖位于“elev”中坐标为 0.1、0.1 的点的单元格的值。 返回的对象是一个生成器：

src_elev.sample([(0.1, 0.1)])

<generator object sample_gen at 0x000001AA0BD6D890>

如果我们想要一次获得所有值，我们可以应用“list”。 结果是“16”：

list(src_elev.sample([(0.1, 0.1)]))

[array([16], dtype=uint8)]

栅格对象也可以是另一个栅格对象的子集，如下面的代码块所示：

…

# ...

空间子集的另一个常见用例是使用布尔掩码，基于具有相同范围和分辨率的另一个栅格或原始栅格，如图所示…为此，我们“擦除”数组中的值 一个栅格的，根据另一个相应的“掩模”栅格。 例如，让我们读取“elev.tif”栅格数组：

elev = src_elev.read(1)
elev

array([[ 1,  2,  3,  4,  5,  6],
       [ 7,  8,  9, 10, 11, 12],
       [13, 14, 15, 16, 17, 18],
       [19, 20, 21, 22, 23, 24],
       [25, 26, 27, 28, 29, 30],
       [31, 32, 33, 34, 35, 36]], dtype=uint8)

并创建相应的random mask：

np.random.seed(1)
mask = np.random.choice([True, False], elev.shape)
mask

array([[False, False,  True,  True, False, False],
       [False, False, False,  True,  True, False],
       [ True, False, False,  True,  True, False],
       [ True,  True,  True, False,  True,  True],
       [False,  True,  True,  True, False,  True],
       [ True,  True, False, False, False, False]])

在上面的代码块中，我们创建了一个名为“mask”的掩码对象，其值随机分配给“True”和“False”。 接下来，我们希望将“elev”的那些值为“False”的值保留在“mask”中（即，它们“不”被屏蔽）。 换句话说，我们想用“mask”屏蔽“elev”。 结果存储在名为“elev1”的副本中。 为了能够在栅格中存储“np.nan”，我们还需要将其转换为浮点型（请参阅@sec-summarizing-raster-objects）：

elev1 = elev.copy()
elev1 = elev1.astype('float64')
elev1[mask] = np.nan
elev1

array([[ 1.,  2., nan, nan,  5.,  6.],
       [ 7.,  8.,  9., nan, nan, 12.],
       [nan, 14., 15., nan, nan, 18.],
       [nan, nan, nan, 22., nan, nan],
       [25., nan, nan, nan, 29., nan],
       [nan, nan, 33., 34., 35., 36.]])

结果显示在@fig-raster-subset 中。

#| label: fig-raster-subset
#| fig-cap: Original raster (left). Raster mask (middle). Output of masking a raster (right).

fig, axes = plt.subplots(ncols=3, figsize=(9,5))
show(elev, ax=axes[0])
show(mask, ax=axes[1])
show(elev1, ax=axes[2])
axes[0].set_title("Original")
axes[1].set_title("Mask")
axes[2].set_title("Result");

在这里插入图片描述

上述方法也可用于用“nan”替换某些值（例如，预期错误的值）：

elev1 = elev.copy()
elev1 = elev1.astype('float64')
elev1[elev1 < 20] = np.nan
elev1

array([[nan, nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan, nan],
       [nan, 20., 21., 22., 23., 24.],
       [25., 26., 27., 28., 29., 30.],
       [31., 32., 33., 34., 35., 36.]])

这些运算实际上是布尔本地运算，因为我们按像元比较两个栅格。 下一小节将更详细地探讨这些操作和相关操作。

4.2 地图代数

“地图代数”这个术语是在 20 世纪 70 年代末创造的，用来描述用于分析地理栅格和（尽管不太突出）矢量数据的“一组约定、功能和技术”（Tomlin 1994）。 在这种情况下，我们将地图代数更狭义地定义为修改或汇总栅格像元值的运算，并参考周围像元、区域或适用于每个像元的统计函数。

地图代数运算往往很快，因为栅格数据集仅隐式存储坐标，因此有句老话“栅格更快，但矢量是校正器”。 栅格数据集中像元的位置可以通过使用其矩阵位置以及数据集的分辨率和原点（存储在标头中）来计算。 然而，对于处理来说，只要我们确保处理后单元格位置仍然相同，单元格的地理位置就几乎不相关。 此外，如果两个或多个栅格数据集共享相同的范围、投影和分辨率，则可以将它们视为矩阵进行处理。

这就是地图代数与 terra 包配合使用的方式。 首先，查询并检查栅格数据集的标头（如果地图代数运算适用于多个数据集），以确保数据集兼容。 其次，地图代数保留了所谓的一对一位置对应关系，这意味着细胞不能移动。 这与矩阵代数不同，在矩阵代数中，值会改变位置，例如在矩阵相乘或相除时。

地图代数（或使用栅格数据的制图建模）将栅格操作分为四个子类（Tomlin 1990），每个子类同时在一个或多个网格上工作：

  • 本地或每单元操作
  • 焦点或邻域操作。 大多数情况下，输出单元值是 3 x 3 输入单元块的结果
  • 区域操作与焦点操作类似，但计算新值的周围像素网格可以具有不规则的大小和形状
  • 全局或逐栅格操作； 这意味着输出像元可能从一个或多个完整栅格中获取其值

这种类型学根据每个像素处理步骤使用的单元数量和输出类型对地图代数运算进行分类。 为了完整起见，我们应该提到，栅格操作也可以按地形、水文分析或图像分类等学科进行分类。 以下部分参考工作示例解释了如何使用每种类型的地图代数运算。

4.3 本地操作

局部操作包括一层或多层中的所有逐单元操作。
栅格代数是局部运算的经典用例 - 这包括从栅格中添加或减去值、对栅格进行平方和乘法。
栅格代数还允许逻辑运算，例如查找大于特定值（在下面的示例中为 5）的所有栅格单元。
使用“numpy”数组操作语法应用本地操作，如下所示：

首先，我们需要读取栅格值：

elev = src_elev.read(1)
elev

array([[ 1,  2,  3,  4,  5,  6],
       [ 7,  8,  9, 10, 11, 12],
       [13, 14, 15, 16, 17, 18],
       [19, 20, 21, 22, 23, 24],
       [25, 26, 27, 28, 29, 30],
       [31, 32, 33, 34, 35, 36]], dtype=uint8)

现在，可以应用任何按元素的数组运算。 例如：

elev + elev

array([[ 2,  4,  6,  8, 10, 12],
       [14, 16, 18, 20, 22, 24],
       [26, 28, 30, 32, 34, 36],
       [38, 40, 42, 44, 46, 48],
       [50, 52, 54, 56, 58, 60],
       [62, 64, 66, 68, 70, 72]], dtype=uint8)

@fig-raster-local-operations 演示了更多示例。

#| label: fig-raster-local-operations
#| fig-cap: 'Examples of different local operations of the elev raster object: adding two rasters, squaring, applying logarithmic transformation, and performing a logical operation.'

fig, axes = plt.subplots(ncols=4, figsize=(9,5))
show(elev + elev, ax=axes[0], cmap="Oranges")
show(elev ** 2, ax=axes[1], cmap="Oranges")
show(np.log(elev), ax=axes[2], cmap="Oranges")
show(elev > 5, ax=axes[3], cmap="Oranges")
axes[0].set_title("elev+elev")
axes[1].set_title("elev ** 2")
axes[2].set_title("np.log(elev)")
axes[3].set_title("elev > 5");

在这里插入图片描述

本地操作的另一个很好的例子是将数值区间分类为组，例如将数字高程模型分为低海拔（第 1 类）、中海拔（第 2 类）和高海拔（第 3 类）。 在这里，我们分配 0-12、12-24 和 24-36 范围内的栅格值，并分别重新分类为值 1、2 和 3。

recl = elev.copy()
recl[(elev > 0)  & (elev <= 12)] = 1
recl[(elev > 12) & (elev <= 24)] = 2
recl[(elev > 24) & (elev <= 36)] = 3

重新分类的结果显示在@fig-raster-reclassify中。

#| label: fig-raster-reclassify
#| fig-cap: Reclassifying a continuous raster into three categories.

fig, axes = plt.subplots(ncols=2, figsize=(9,5))
show(elev, ax=axes[0], cmap="Oranges")
show(recl, ax=axes[1], cmap="Oranges")
axes[0].set_title("Original")
axes[1].set_title("Reclassified");

在这里插入图片描述

归一化植被指数（NDVI）的计算是众所周知的局部（逐像素）栅格操作。 它返回一个值介于 -1 和 1 之间的栅格； 正值表示存在活植物（大多数> 0.2）。 NDVI 是根据遥感图像的红色和近红外 (NIR) 波段计算得出的，通常来自 Landsat 或 Sentinel 等卫星系统。 植被大量吸收可见光谱中的光，尤其是红色通道中的光，同时反射近红外光，解释了 NVDI 公式：

N D V I = N I R − R e d N I R + R e d NDVI=\frac{NIR-Red} {NIR+Red} NDVI=NIR+RedNIR−Red​

让我们计算锡安国家公园多光谱卫星文件的 NDVI。

multi_rast = src_multi_rast.read()
nir = multi_rast[3,:,:]
red = multi_rast[2,:,:]
ndvi = (nir-red)/(nir+red)

将值 >1 转换为“无数据”：

ndvi[ndvi>1] = np.nan

当使用“show”函数绘制 RGB 图像时，该函数假设：

  • 浮点数的值在“[0,1]”范围内，整数的值在“[0,255]”范围内（否则会被截断）
  • 波段顺序为 RGB

为了“准备”用于“show”的多波段栅格，我们反转了波段的顺序（最初是 BGR+NIR），并除以最大值以将最大值设置为“1”：

multi_rast_rgb = multi_rast[(2,1,0), :, :] / multi_rast.max()

结果显示在@fig-raster-ndvi 中。

#| label: fig-raster-ndvi
#| fig-cap: RGB image (left) and NDVI values (right) calculated for the example satellite file of the Zion National Park.

fig, axes = plt.subplots(ncols=2, figsize=(9,5))
show(multi_rast_rgb, ax=axes[0], cmap="RdYlGn")
show(ndvi, ax=axes[1], cmap="Greens")
axes[0].set_title("RGB image")
axes[1].set_title("NDVI");

在这里插入图片描述

4.4 焦点操作

虽然局部函数在一个单元上运行，但可能来自多个层，但焦点操作会考虑中心（焦点）单元及其邻居。
所考虑的邻域（也称为内核、过滤器或移动窗口）通常大小为 3×3 个单元（即中心单元及其八个周围的邻居），但可以采用任何其他（不一定是矩形）形状，如 由用户定义。
焦点操作将聚合函数应用于指定邻域内的所有单元，使用相应的输出作为中心单元的新值，然后移动到下一个中心单元（图…）。
此操作的其他名称是空间过滤和卷积（Burrough、McDonnell 和 Lloyd 2015）。

在Python中，scipy.ndimage包有一个全面的函数集合来执行numpy数组的过滤， 例如：

  • 最小过滤器
  • 最大过滤器
  • uniform_filter （即均值过滤器）
  • 中值过滤器
    ETC。

在这组函数中，我们使用以下之一定义移动窗口的形状：

  • size——单个数字或元组，意味着这些维度的过滤器
  • footprint——一个布尔数组，代表窗口形状和所包含元素的标识

除了特定的内置过滤器外，

  • convolve 在乘以自定义的 weights 数组后应用求和函数
  • generic_filter 可以传递任何自定义函数，用户可以在其中指定任何类型的基于窗口的自定义计算。

例如，这里我们在“elev”上应用窗口大小为“3”的最小过滤器：

elev

array([[ 1,  2,  3,  4,  5,  6],
       [ 7,  8,  9, 10, 11, 12],
       [13, 14, 15, 16, 17, 18],
       [19, 20, 21, 22, 23, 24],
       [25, 26, 27, 28, 29, 30],
       [31, 32, 33, 34, 35, 36]], dtype=uint8)

elev_min = scipy.ndimage.minimum_filter(elev, size=3)
elev_min

array([[ 1,  1,  2,  3,  4,  5],
       [ 1,  1,  2,  3,  4,  5],
       [ 7,  7,  8,  9, 10, 11],
       [13, 13, 14, 15, 16, 17],
       [19, 19, 20, 21, 22, 23],
       [25, 25, 26, 27, 28, 29]], dtype=uint8)

应特别注意边缘像素。 应该如何计算？ scipy.ndimage 通过 mode 参数提供了几个选项：

  • reflect（默认）
  • 常数
  • 最近的
    *镜子
  • 包裹

有时人为地延伸光栅边缘被认为是不合适的。
换句话说，我们可能希望生成的栅格仅包含具有“完整”窗口的像素值，例如具有统一的样本大小或因为所有方向上的值都很重要（例如在地形计算中）。
在scipy.ndimage中没有特定的选项不来扩展边缘。
但是，为了获得相同的效果，可以根据过滤器大小在多个行和列中为过滤后的数组的边缘分配“nan”。
例如，当使用“size=3”的过滤器时，像素的第一“层”可能会被分配为“nan”，反映这些像素具有不完整的 3*3 邻域的事实：

elev_min = elev_min.astype('float')
elev_min[:, [0, -1]] = np.nan
elev_min[[0, -1], :] = np.nan
elev_min

array([[nan, nan, nan, nan, nan, nan],
       [nan,  1.,  2.,  3.,  4., nan],
       [nan,  7.,  8.,  9., 10., nan],
       [nan, 13., 14., 15., 16., nan],
       [nan, 19., 20., 21., 22., nan],
       [nan, nan, nan, nan, nan, nan]])

我们可以快速检查输出是否符合我们的预期。
在我们的示例中，最小值必须始终是移动窗口的左上角（请记住，我们是通过从左上角开始按行将像元值递增 1 创建输入栅格的）。

焦点函数或滤波器在图像处理中起着主导作用。 低通或平滑滤波器使用均值函数来消除极值。
对于分类数据，我们可以用众数替换均值，众数是最常见的值。
相比之下，高通滤波器强调特征。 此处以线检测拉普拉斯滤波器和索贝尔滤波器为例。

地形处理，即坡度、坡向和流向等地形特征的计算，依赖于焦点函数。
richdem 包中的 TerrainAttribute 函数可用于计算常见的指标，通过 attrib 参数，即：

  • slope_riserun Horn (1981) doi: 10.1109/PROC.1981.11918
  • slope_percentage Horn (1981) doi: 10.1109/PROC.1981.11918
  • slope_degrees Horn (1981) doi: 10.1109/PROC.1981.11918
  • slope_radians Horn (1981) doi: 10.1109/PROC.1981.11918
  • aspect Horn (1981) doi: 10.1109/PROC.1981.11918
  • curvature Zevenbergen and Thorne (1987) doi: 10.1002/esp.3290120107
  • planform_curvature Zevenbergen and Thorne (1987) doi: 10.1002/esp.3290120107
  • profile_curvature Zevenbergen and Thorne (1987) doi: 10.1002/esp.3290120107

4.5 Zonal operations

就像焦点操作一样，区域操作将聚合函数应用于多个栅格单元。
然而，第二个栅格通常具有分类值，在区域操作的情况下定义区域过滤器（或“区域”），这与上一节中介绍的焦点操作情况下的预定义邻域窗口相反。
因此，定义区域过滤器的栅格单元不一定必须是邻居。
我们的粒度栅格就是一个很好的例子，如图 3.2 的右图所示：不同的粒度在整个栅格中不规则地分布。
最后，区域操作的结果是按区域分组的汇总表，这就是为什么该操作在 GIS 领域也称为区域统计。 这与返回光栅对象的焦点操作形成对比。

为了进行演示，让我们回到“grain”和“elev”栅格（图 3.2）。 为了计算区域统计数据，我们使用具有栅格值的数组。 elev 数组之前已经导入：

grain = src_grain.read(1)
grain

array([[1, 0, 1, 2, 2, 2],
       [0, 2, 0, 0, 2, 1],
       [0, 2, 2, 0, 0, 2],
       [0, 0, 1, 1, 1, 1],
       [1, 1, 1, 2, 1, 1],
       [2, 1, 2, 2, 0, 2]], dtype=uint8)

我们的目的是计算由 grain 值定义的每个区域中 elevation 的平均值（或任何其他汇总函数）。 首先，我们可以使用“np.unique”获取定义区域的唯一值：

np.unique(grain)

array([0, 1, 2], dtype=uint8)

现在，我们可以使用dictionary conprehension将“elev”数组“拆分”为单独的一维数组，每个“grain”组都有值 ，键是唯一的“grain”值：

z = {i: elev[grain == i] for i in np.unique(grain)}
z

{0: array([ 2,  7,  9, 10, 13, 16, 17, 19, 20, 35], dtype=uint8),
 1: array([ 1,  3, 12, 21, 22, 23, 24, 25, 26, 27, 29, 30, 32], dtype=uint8),
 2: array([ 4,  5,  6,  8, 11, 14, 15, 18, 28, 31, 33, 34, 36], dtype=uint8)}

在此阶段，我们可以扩展字典理解表达式来计算与每个粒度类别相关的平均海拔。 我们不是将高程值 (elev[grain==i]) 放入字典值中，而是将它们的平均值 (elev[grain==i].mean()) 放入：

z = {i: elev[grain == i].mean() for i in np.unique(grain)}
z

{0: 14.8, 1: 21.153846153846153, 2: 18.692307692307693}

这将返回每个类别的统计数据，这里是每个粒度类别的平均海拔。 例如，以粒度“0”为特征的像素的平均高程为“14.8”，依此类推。

4.6 合并栅格

假设我们想要计算 NDVI（请参阅@sec-raster-local-operations），并且还想要根据研究区域内观测的高程数据计算地形属性。
此类计算依赖于遥感信息。
相应的图像通常被划分为覆盖特定空间范围的场景（即“图块”），并且通常一个研究区域覆盖多个场景。
然后，我们需要合并（也称为“马赛克”）我们研究区域覆盖的场景。
在最简单的情况下，我们可以合并这些场景，即将它们并排放置。
然而，无论场景数量多少，步骤都是相同的。
如果所有场景都“对齐”（即共享相同的原点和分辨率），则可以将其视为简单地将它们粘合到一个大光栅中； 否则，所有场景都将重新采样（请参阅@sec-raster-resampling）到第一个场景定义的网格。

例如，让我们合并来自奥地利 (aut.tif) 和瑞士 (ch.tif) 两个 SRTM 高程图块的数字高程数据。
可以使用函数“rasterio.merge.merge”完成合并，该函数接受光栅文件连接的“列表”，并返回新的“ndarray”和“转换”：

src_1 = rasterio.open('data/aut.tif')
src_2 = rasterio.open('data/ch.tif')
out_image, out_transform = rasterio.merge.merge([src_1, src_2])

输入和结果都显示在 @fig-raster-merge 中：

#| label: fig-raster-merge
#| fig-cap: Raster merging

fig, axes = plt.subplots(2, 2, figsize=(8,4))
show(src_1, ax=axes[0][0])
show(src_2, ax=axes[0][1])
show(out_image, transform=out_transform, ax=axes[1][0])
axes[1][1].axis('off')
axes[0][0].set_title('aut.tif')
axes[0][1].set_title('ch.tif')
axes[1][0].set_title('Mosaic (aut.tif+ch.tif)');

在这里插入图片描述

默认情况下（method='first'），重叠区域保留first 栅格的值。 其他可能的方法有：

  • 'last'—最后一个栅格的值
  • 'min'—最小值
  • 'max'—最大值

当处理非重叠图块时，例如“aut.tif”和“ch.tif”（上面），“method”参数没有实际效果。 然而，当我们想要组合不同日期拍摄的场景的光谱图像时，它就变得相关了。 上述四个“方法”选项并没有涵盖我们想要计算平均值时通常需要的场景，例如从一组部分重叠的卫星图像（例如 Landsat）计算季节性平均 NDVI 图像。 用于计算马赛克以及“平均”任何重叠的“rasterio.merge.merge”的替代工作流程可能需要执行两个步骤：

  • 将所有场景重新采样到一个通用的“全局”网格中（@sec-raster-resampling），从而生成一系列“匹配”栅格（每个场景周围的区域设置为“无数据”）
  • 通过栅格代数 (@sec-raster-local-operations) 对栅格进行平均，如 np.mean(m,axis=0) 或 np.nanmean(m,axis=0)，其中 m 是多波段数组，它将返回单波段平均值数组
    2]], dtype=uint8)

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/675509.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

文件搜索引擎的搭建Elasticsearch+Fscrawler+SearchUI+Git+Nginx

文件搜索引擎的搭建Elasticsearch+Fscrawler+SearchUI+Git+Nginx

文章目录 前言如何搭建文档搜索引擎服务器架构环境准备一、搭建Elasticsearch二、搭建Fscrawler三、搭建SearchUI服务四、定时拉取Git文件五、搭建Nginx文件下载服务器 前言 搭建一套文档搜索引擎。有时候&#xff0c;我们有一批文档&#xff0c;需要在这批文档中查找想要的内…
阅读更多...
【从零开始学习JAVA | 第十篇】StringBuild介绍

【从零开始学习JAVA | 第十篇】StringBuild介绍

前言&#xff1a; 本文会对StringBuild类进行详细的介绍&#xff0c;他相比较于String&#xff0c;在进行字符串拼接的时候大大提高了效率&#xff0c;是一个 很实用的工具。 StringBulider&#xff1a; StringBuilder是Java中的一个类&#xff0c;用于在一次创建和初始化后&a…
阅读更多...
Java中反射机制,枚举,Lambda的使用

Java中反射机制,枚举,Lambda的使用

目录 一、反射机制 1、含义 2、作用 3、※反射相关的几个类 3.1、Class类&#xff08;Class对象是反射的基石&#xff09; 3.2、Class类中相关的方法 3.2.1 (※重要)常用获得类相关的方法 3.2.2 (※重要)常用获得类中属性、变量Field相关的方法 3.2.3 获得类中注解相…
阅读更多...
详细设计报告

详细设计报告

聊天系统设计与实现详细设计报告 1.编写目的 详细设计的主要任务是概要设计方案做完善和细化，本阶段主要对聊天应用系统进行过程化的描述，详细确定每一个功能模块的实现方式、执行流程，为程序员编码提供依据。设计用户界面。 2.总体方案确认 (1)系统总体结构确认 该项目管理…
阅读更多...
Redis原理 - 通信协议RESP

Redis原理 - 通信协议RESP

原文首更地址&#xff0c;阅读效果更佳&#xff01; Redis原理 - 通信协议RESP | CoderMast编程桅杆https://www.codermast.com/database/redis/redis-communication-protocol.html RESP协议 Redis 是一个 CS 架构的软件&#xff0c;通信一般分两步&#xff08;不包括pipeli…
阅读更多...
计算机组成与设计Patterson Hennessy 笔记_1 计算机概要与技术

计算机组成与设计Patterson Hennessy 笔记_1 计算机概要与技术

Patterson & Hennessy 计算机概要与技术 计算机应用包括&#xff1a;个人计算机PC&#xff0c;服务器&#xff0c;嵌入式计算机。后PC时代出现了个人移动设备PMD&#xff08;手机&#xff09;&#xff0c;云计算&#xff08;在网络上提供服务的大服务器集群&#xff0c;供…
阅读更多...
Unity Mac最新打苹果包流程

Unity Mac最新打苹果包流程

作者介绍&#xff1a;铸梦xy。IT公司技术合伙人&#xff0c;IT高级讲师&#xff0c;资深Unity架构师&#xff0c;铸梦之路系列课程创始人。 IOS详细打包流程1.申请APPID2.申请开发证书3.创建描述文件 IOS详细打包流程 1.申请AppID 2.创建证书 3.申请配置文件&#xff08;又名描…
阅读更多...
NodeJS KOA⑩②

NodeJS KOA⑩②

文章目录 ✨文章有误请指正&#xff0c;如果觉得对你有用&#xff0c;请点三连一波&#xff0c;蟹蟹支持&#x1f618;前言KOA Koa vs Express Koa更轻量 Koa~Context对象 Koa~异步流程控制 Koa~中间件模型Koa路由 1.1基本使用 2.2请求方式2.2.1规范写法2…
阅读更多...
数据结构——带头节点的双向循环列表

数据结构——带头节点的双向循环列表

带头节点的双向循环链表是一种特殊的双向链表&#xff0c;它与普通的双向链表相比&#xff0c;最大的区别是链表头结点的 next 指针不再指向第一个实际节点&#xff0c;而是指向链表中的第一个节点。同时&#xff0c;链表尾结点的 prev 指针也不再指向 NULL&#xff0c;而是指向…
阅读更多...
轻松配置深度学习模型 ?

轻松配置深度学习模型 ?

动动发财的小手&#xff0c;点个赞吧&#xff01; 由于所有模块都需要大量参数和设置&#xff0c;因此管理深度学习模型可能很困难。训练模块可能需要诸如 batch_size 或 num_epochs 之类的参数或学习率调度程序的参数。同样&#xff0c;数据预处理模块可能需要 train_test_spl…
阅读更多...
Java接口幂等性，如何重试？

Java接口幂等性，如何重试？

Java接口幂等性&#xff0c;如何重试&#xff1f; 文章目录 Java接口幂等性&#xff0c;如何重试&#xff1f;前言一、幂等性是什么&#xff1f;二、为什么要幂等性&#xff1f;三、使用什么办法实现幂等性&#xff1f;1.insert前先select2.加悲观锁3.加乐观锁4.加唯一索引5.Re…
阅读更多...
uniapp系列-uni.getAppBaseInfo() versionCode appVersion 值不对应该怎么解决？

uniapp系列-uni.getAppBaseInfo() versionCode appVersion 值不对应该怎么解决？

今天看到一个BUG 问题描述 我们使用uniapp的官方文档中uni.getAppBaseInfo()后获取的 appVersionCode appVersion &#xff0c;发现获得的结果和我们实际设置的不一致&#xff0c;不是manifest.json里面的值&#xff0c;如下图所示官方文档&#xff1a;https://uniapp.dcloud…
阅读更多...
还原大师（MD5）

还原大师（MD5）

根据题目提示&#xff0c;都猜得到这应该跟MD5的加密形式有关系 我好像还没有具体了解过MD5编码的格式&#xff0c;或许本题可以通过MD5的编码格式推导出字符串 但是说实话&#xff0c;MD5的加密方式没有找到详细简介的文章 然后我就去网上百度了一下&#xff0c;经过大佬wp的洗…
阅读更多...
java springboot整合MyBatis演示增删查改操作

java springboot整合MyBatis演示增删查改操作

前面我的文章 java springboot整合MyBatis做数据库查询操作讲述了整合springboot整合MyBatis 做了根据id查询的语句 那么 我们现在按它搭建的项目继续 我们在staffDao中添加一个insert函数 参考代码如下 Insert("insert into staff(name, age, status, departmentid) va…
阅读更多...
chatgpt赋能python：Python编程实现1+22+333，解密方法

chatgpt赋能python：Python编程实现1+22+333，解密方法

Python编程实现122333&#xff0c;解密方法 在Python编程开发中&#xff0c;我们经常需要求解不同类型的算数表达式&#xff0c;其中求解一系列类似122333的表达式是一个比较常见的需求。本文将会介绍如何使用Python语言方便地求解这类表达式&#xff0c;为大家提供一种针对此…
阅读更多...
Redis中AOF和RDB

Redis中AOF和RDB

在Redis的持久化中&#xff0c;常使用的两个手段便是AOF和RDB进行持久化。 RDB&#xff08;Redis DataBase&#xff09;是Redis的持久化方式之一&#xff0c;在配置文件中&#xff0c;我们可以找到 对Redis进行持久化配置&#xff0c;而RDB在持久化时是怎么样进行工作的呢&…
阅读更多...
ARM、ARM架构、ARM架构芯片

ARM、ARM架构、ARM架构芯片

ARM是一种基于精简指令集&#xff08;RISC&#xff09;的处理器架构&#xff0c;它由英国的ARM公司设计和授权。 ARM芯片具有低功耗、高性能、高集成度等特点&#xff0c;广泛应用于嵌入式系统、移动设备、物联网、服务器等领域。本文将介绍ARM的各类芯片&#xff0c;包括其特…
阅读更多...
Java-API简析_java.lang.ClassLoader类（基于 Latest JDK）（浅析源码）

Java-API简析_java.lang.ClassLoader类（基于 Latest JDK）（浅析源码）

【版权声明】未经博主同意&#xff0c;谢绝转载&#xff01;&#xff08;请尊重原创&#xff0c;博主保留追究权&#xff09; https://blog.csdn.net/m0_69908381/article/details/131345825 出自【进步*于辰的博客】 其实我的【Java-API】专栏内的博文对大家来说意义是不大的。…
阅读更多...
CentOS 7.9 安装 Jenkins

CentOS 7.9 安装 Jenkins

CentOS 7.9 安装 Jenkins 文章目录 CentOS 7.9 安装 Jenkins一、概述二、安装1、安装 OpenJDK2、安装 Jenkins3、启动 Jenkins4、给 Jenkins 放行端口 三、初始化 Jenkins 配置1、访问2、解锁 Jenkins3、配置清华大学的源地址4、安装插件5、创建管理员用户6、完成安装 四、功能…
阅读更多...
TypeScript ~ TS 掌握自动编译命令 ③

TypeScript ~ TS 掌握自动编译命令 ③

作者 : SYFStrive 博客首页 : HomePage &#x1f4dc;&#xff1a; TypeScript ~ TS &#x1f4cc;&#xff1a;个人社区&#xff08;欢迎大佬们加入&#xff09; &#x1f449;&#xff1a;社区链接&#x1f517; &#x1f4cc;&#xff1a;觉得文章不错可以点点关注 &…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023