使用-GeoPandas-实现地图边界提取
使用 GeoPandas 实现地图边界提取
原文：
towardsdatascience.com/easy-map-boundary-extraction-with-geopandas-7ee1465e1344/
快速成功数据科学
提取 DALL-E3 想象的边界
我最近遇到了一个地理游戏，游戏中展示了国家之间的一段边界，并要求你命名它们。这里有一个例子；你能猜出共享这条边界的两个国家吗？
猜猜这个著名的 "开放" 边界（作者）
这里还有一个 "开放" 边界。北边是哪个国家？
猜猜这个边界（作者）
哪些国家共享这条超级长的边界：
猜猜这个边界！（作者）
老一代的婴儿潮人群可能熟悉这个臭名昭著的边界地区：
猜猜这个边界（作者）
按顺序，这些是美墨、美加、阿根廷智利和越南老挝柬埔寨边界。
除了测试你的地理知识外，你还可以使用提取的边界制作有趣的信息图表。在这里，我展示了四个最长共享边界的长度：
作者认为最长的 4 个共享边界
在这个快速成功的数据科学项目中，我们将使用 Python 的开源 GeoPandas 地理空间库来提取、绘制和计算政治边界的长度。这个过程也可以用于存储在 shapefile 中的其他边界，例如省、州、县和市界。
流程
这里是我们将使用的高级流程：
-
加载 shapefile： 使用 GeoPandas 加载包含国家边界的 shapefile。shapefile 是地理信息系统（GIS）软件的地理空间矢量数据格式。
-
筛选相关国家： 使用 GeoPandas 将感兴趣的国家的边界筛选到单独的 GeoPandas GeoDataFrames 中。
-
提取国家边界： 在 GeoDataFrames 的 "Geometry" 列中定位国家边界。此列包含用于绘制边界的坐标。
-
查找共享边界： 通过计算主要国家边界与相邻国家边界的交集来识别共享边界。
-
创建边界坐标列表： 将边界数据从几何对象（例如，
LineString或MultiLineString）转换为列表。 -
绘制边界： 将边界作为 GeoPandas 的
GeoSeries绘制。
我们还将探讨如何计算边界的长度。
有关 shapefiles 和 GeoPandas 的概述，包括安装说明，请参阅此文章：
下载 Shapefile
Natural Earth是一个公共领域数据集，包括 shapefiles。您可以从这里下载国家边界shapefile（见下图）。这是 Natural Earth 在 1:10,000,000（1 厘米=100 公里）比例尺上最详细的数据集。
点击黄色高亮的下载链接（Natural Earth 网页）
上一张图中高亮的链接将下载 _ne_10m_admin_0countries.zip文件。没有必要解压它，因为 GeoPandas 可以处理压缩数据。为了方便，将其存储在 Python 脚本或笔记本的同一文件夹中。
完整代码
完整的注释代码如下。我将在接下来的章节中更详细地描述它。
import matplotlib.pyplot as plt
import geopandas as gpd
# Define the primary country name:
# (Must match names in shapefile).
COUNTRY1_NAME = 'Venezuela'
# Define the name(s) of adjacent countries:
# (Must match names in shapefile).
COUNTRY2_NAME = ['Colombia', 'Brazil', 'Guyana']
# Load the shapefile as a GeoDataFrame:
gdf = gpd.read_file('ne_10m_admin_0_countries.zip')
gdf = gdf.rename(columns={'ADMIN': 'name'})
# Make separate GDF for country1 and adjacent countries:
country1 = gdf[gdf['name'] == COUNTRY1_NAME]
others = gdf[gdf['name'].isin(COUNTRY2_NAME)]
# Extract the boundaries:
country1_border = country1.boundary.unary_union
others_border = others.boundary.unary_union
# Find intersections
shared_borders = country1_border.intersection(others_border)
# Ensure geometry is a MultiLineString or LineString:
if shared_borders.geom_type == 'GeometryCollection':
borders = [geom for geom in shared_borders.geoms if
geom.geom_type in ['LineString', 'MultiLineString']]
elif shared_borders.geom_type in ['MultiLineString', 'LineString']:
borders = [shared_borders]
else:
borders = []
# Create a dynamic title:
countries_in_title = f"{COUNTRY1_NAME} and {', '.join(COUNTRY2_NAME)}"
title = f"Shared Borders Between {countries_in_title}"
# Set up the figure:
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_title(title, fontsize=14)
# ax.set_xticks([])
# ax.set_yticks([])
# ax.set_xticklabels([])
# ax.set_yticklabels([])
# Plot the borders:
gpd.GeoSeries(borders).plot(ax=ax,
color='firebrick',
linewidth=2)
# plt.savefig(f"{COUNTRY1_NAME}_border_segment.png", dpi=600)
plt.show()
此代码绘制了委内瑞拉、哥伦比亚和圭亚那的共享边界：
委内瑞拉、哥伦比亚和圭亚那之间的边界（作者提供）
导入库和分配常量
首先，导入 Matplotlib（用于绘图）和 GeoPandas（用于处理 shapefiles）。接下来，我们为主要的国名字符串（COUNTRY1_NAME）和共享边界的相邻国家名称（COUNTRY2_NAME）分配常量。后者变量使用一个包含一个或多个国家名称的 Python 列表。
import matplotlib.pyplot as plt
import geopandas as gpd
# Define the primary country name:
# (Must match names in shapefile).
COUNTRY1_NAME = 'Venezuela'
# Define the name(s) of adjacent countries:
# (Must match names in shapefile).
COUNTRY2_NAME = ['Colombia', 'Brazil', 'Guyana']
注意：用作常量的国家名称必须与 shapefile 中的名称匹配。我们将在下一节中更详细地讨论这个问题。
加载和准备 Shapefiles
现在我们将下载的 shapefile 作为名为"gdf"的 GeoDataFrame 加载。注意我们不需要解压 shapefile。
# Load the shapefile as a GeoDataFrame:
gdf = gpd.read_file('ne_10m_admin_0_countries.zip')
gdf = gdf.rename(columns={'ADMIN': 'name'})
# Make separate GDF for country1 and adjacent countries:
country1 = gdf[gdf['name'] == COUNTRY1_NAME]
others = gdf[gdf['name'].isin(COUNTRY2_NAME)]
Natural Earth 的 shapefile 使用无意义的名称"ADMIN"作为国家名称列。我们可以使用 GeoPandas 的rename()方法来修复这个问题。
接下来，我们使用新创建的"名称"列创建两个新的 GeoDataFrames。对于第一个（country1），我们使用 COUNTRY1_NAME 进行筛选。由于 COUNTRY2_NAME 变量是一个列表，我们添加了一个转折，并使用了isin()方法。这样，列表中任何匹配的名称都将用于筛选。
说到匹配的名称，如果您不确定国家名称的拼写或存储在 GeoDataFrame 中（例如，“United States of America”与“USA”），您可以使用此命令打印所有名称：
print(gdf['name'].unique())
寻找共享边界
在 GeoDataFrames 中加载数据后，我们现在提取每个国家的边界并找到它们的交点。
# Extract the boundaries:
country1_border = country1.boundary.unary_union
others_border = others.boundary.unary_union
# Find intersections
shared_borders = country1_border.intersection(others_border)
# Create list from geometry object:
if shared_borders.geom_type == 'GeometryCollection':
borders = [geom for geom in shared_borders.geoms if
geom.geom_type in ['LineString', 'MultiLineString']]
elif shared_borders.geom_type in ['MultiLineString', 'LineString']:
borders = [shared_borders]
else:
borders = []
要提取边界，我们将 GeoPandas 的boundary操作与 Shapely 的unary_union属性（Shapely 是一个库，GeoPandas 在底层用于几何操作）链式使用。第一部分产生一个包含 Shapely LineString或MultiLineString对象的GeoSeries。第二部分通过在GeoSeries中的所有几何体上执行并集操作来聚合多个几何体（例如，从country1.boundary）。它返回一个单一的 Shapely Geometry对象（例如，LineString，MultiLineString或GeometryCollection）。
下一步使用 Shapely 的intersection()方法来识别共享边界。此方法计算两个几何体重叠的部分，例如线或多边形。它根据输入几何体及其空间关系返回许多 Shapely 几何体对象之一。例如，线字符串、点、多边形等。
由于intersection()方法可能输出的结果复杂，我们必须检查数据类型是否有效，以便与线进行工作。我们通过一个条件块来完成这项工作，该块返回最终边界坐标作为 Python 列表中的 Shapely 几何体。
绘制结果
代码的其余部分用于绘制边界。在设置完图形后，我们将borders从列表转换为 GeoSeries（gpd.GeoSeries(borders)），然后在这个对象上调用plot()。
# Create a dynamic title:
countries_in_title = f"{COUNTRY1_NAME} and {', '.join(COUNTRY2_NAME)}"
title = f"Shared Borders Between {countries_in_title}"
# Set up the figure:
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_title(title, fontsize=14)
# Uncomment to hide axis labels and ticks:
# ax.set_xticks([])
# ax.set_yticks([])
# ax.set_xticklabels([])
# ax.set_yticklabels([])
# Plot the borders:
gpd.GeoSeries(borders).plot(ax=ax,
color='firebrick',
linewidth=2)
# plt.savefig(f"{COUNTRY1_NAME}_border_segment.png", dpi=600)
plt.show()
使用 GeoSeries 进行绘图是首选的，因为 GeoPandas 的绘图功能与地理投影集成，并允许更高级的空间数据处理。
计算边界的长度
理论上，在提取边界后我们可以计算边界的长度。第一步是在加载数据后，将 GeoDataFrame 的坐标参考系统（CRS）转换为等面积投影（例如，EPSG:2163）：
gdf = gdf.to_crs('EPSG:2163')
然后，在将边界作为 Shapely 几何体列表提取出来之后，我们可以运行以下代码：
# Calculate the total length of the shared border:
total_length = sum(line.length for line in borders)
print(f"Total shared border length (in projected CRS units): ")
print(f"{total_length / 1000} km")
borders中的每个几何体（例如，LineString或MultiLineString）都有一个由 Shapely 提供的.length属性。该属性是几何体在当前 CRS 数据单位中的计算长度。
不幸的是，使用我们的 Natural Earth 数据集，我们不会得到准确的结果。根据他们的网站，它被设计为"在大型墙式海报上展示世界"。这意味着它为了简化而平滑了细节。
计算一个国家边界的真实长度是一项棘手的工作。你需要一个非常高分辨率的 shapefile 来处理崎岖的边界。例如，考虑俄罗斯和哈萨克斯坦的边界，它在某些地方沿着蜿蜒的河流：
部分哈萨克斯坦-俄罗斯边界（来源：Google Maps）
大多数形状文件无法捕捉到这条边界的每一个曲折，因此我们会低估其长度。当然，我们永远无法知道任何边界的真实长度，因为我们永远不会有一个具有无限分辨率的形状文件。
除了高分辨率形状文件外，您可能还需要使用多个 CRS。地理 CRS 中的纬度相关畸变使得长度看起来更长或更短，这取决于与极地的距离。跨越超过六度纬度或经度的边界应分成段，并使用一个合适的局部CRS 进行重新投影，以最小化每个段的畸变。选项包括UTM 区域或其他本地优化的投影。
测量方法也同样关键。对于长边界的长度计算应使用大地测量方法，这些方法考虑到了地球的曲率。GIS 工具如 GeoPandas 在使用地理坐标参考系统（CRS）时支持大地测量长度计算。
摘要
GeoPandas 提供了一些工具，可以帮助您找到并提取地图区域之间共享的边界线。一旦提取出来，这些边界线就可以进行绘制和测量。这些测量的准确性取决于形状文件的分辨率以及所使用的投影和测量方法。
感谢！
感谢阅读，请关注我，未来将会有更多关于快速成功数据科学项目的更新。
浙公网安备 33010602011771号