Python GeoPandas基础知识:地图、投影和空间连接

0 阅读14分钟

GeoPandas 扩展了 pandas 的功能,使在 Python 中处理地理空间数据更加直观和强大。如果您希望在 Python 中执行地理空间任务,并需要一个具有类似 pandas API 的库,那么 GeoPandas 是一个绝佳的选择。本教程将向您展示如何完成四个常见的地理空间任务:读取数据、绘制地图、应用投影以及进行空间连接。

看完本教程,您将明白:

  • GeoPandas 扩展了 pandas 的功能,使其支持空间数据。这些数据通常存储在 geometry 列中,并允许进行空间操作,例如投影和空间连接;而Folium则专注于在数据准备后创建更丰富的交互式 Web 地图。
  • 您可以使用 .crs检查 CRS(Coordinate Reference Systems,坐标参考系统) ,并使用 .to_crs() 和权威代码(如 EPSG:4326ESRI:54009重新投影数据。
  • 地理坐标参考系以度为单位存储经度和纬度,而投影坐标参考系使用米或英尺等线性单位进行面积和距离计算。
  • 空间连接使用 .sjoin() 以及谓词,例如 "within""intersects",并且两个输入必须共享相同的 CRS ,否则关系将计算错误。

以下是 GeoPandas 与其他替代库的比较:

Use CasePick pandasPick FoliumPick GeoPandas
Tabular data analysis-
Mapping-
Projections, spatial joins--

GeoPandas 在 pandas 的基础上,增加了对地理空间数据和操作(例如投影空间连接)的支持。它还包含用于创建地图的工具。Folium专注于交互式、基于 Web 的地图,用户可以对其进行更深入的自定义,从而与 GeoPandas 形成互补。

GeoPandas入门

首先,您需要准备环境并加载一个小型数据集,该数据集将在整个教程中使用。在接下来的两个小节中,您将安装必要的软件包并读取纽约市行政区边界的示例数据集。这将为您提供一个具体的 GeoDataFrame,供您在学习核心概念的过程中进行探索。

安装 GeoPandas

本教程使用两个软件包:geopandas 用于处理地理数据,geodatasets 用于加载示例数据。建议将这些软件包安装在虚拟环境中,这样可以确保项目与系统其他部分隔离,并便于管理其依赖项

虚拟环境激活后,您可以install both packages with pip

$ python -m pip install "geopandas[all]" geodatasets

使用 [all] 选项可确保您拥有读取数据、转换坐标系和创建图表所需的一切。对于大多数读者来说,这套方案开箱即用。

如果遇到安装问题,项目维护者在官方安装页面上提供了替代安装选项。

读取和数据

大多数地理空间数据集采用GeoJSONshapefile格式。read_file() 函数可以读取这两种格式,并且接受本地文件路径或 URL 作为参数。

在下面的示例中,您将使用 read_file() 加载纽约市行政区边界 (NYBB) 数据集。geodatasets 包提供了一个便捷的数据集路径,因此您无需手动下载任何内容。您还将删除不必要的列:

>>> import geopandas as gpd
>>> import matplotlib.pyplot as plt
>>> from geodatasets import get_path
>>> path_to_data = get_path("nybb")
>>> nybb = gpd.read_file(path_to_data)
>>> nybb = nybb[["BoroName", "Shape_Area", "geometry"]]
>>> nybb
    BoroName        Shape_Area      geometry
0   Staten Island   1.623820e+09    MULTIPOLYGON (((970217.022 145643.332, ....
1   Queens          3.045213e+09    MULTIPOLYGON (((1029606.077 156073.814, ...
2   Brooklyn        1.937479e+09    MULTIPOLYGON (((1021176.479 151374.797, ...
3   Manhattan       6.364715e+08    MULTIPOLYGON (((981219.056 188655.316, ....
4   Bronx           1.186925e+09    MULTIPOLYGON (((1012821.806 229228.265, ...
>>> type(nybb)
<class 'geopandas.geodataframe.GeoDataFrame'>
>>> type(nybb["geometry"])
<class 'geopandas.geoseries.GeoSeries'>

nybb 是一个GeoDataFrame。GeoDataFrame拥有行、列以及pandas DataFrame的所有方法。区别在于它通常包含一个特殊的 geometry 列,用于存储地理形状,而不是纯数字或文本。

geometry 列是一个 GeoSeries。它的行为类似于普通的 pandas Series,但它的值是空间对象,您可以对其进行映射并运行空间查询。在 nybb 数据集中,每个行政区的几何图形都是一个 MultiPolygon——由多个多边形组成——因为每个行政区都包含多个岛屿。很快,您将使用这些几何图形来创建地图并运行空间操作,例如查找某个点属于哪个行政区。

映射数据

加载 GeoDataFrame 后,最快捷的数据可视化方法之一就是了解数据。在本节中,您将学习如何创建静态地图和交互式地图。这样,您可以检查形状、发现规律,并确认几何图形是否符合预期。

创建静态地图

正如pandas 绘图指南中所述,pandas DataFrame提供了一个 .plot() 方法,用于快速创建散点图或柱状图等可视化图表。GeoDataFrame 类重写了 .plot() 方法,因此,当数据包含几何图形时,结果会生成地图而不是图表。

在这里,您可以调用 nybb.plot() 来绘制纽约市五个行政区的轮廓:

>>> nybb.plot()
<Axes: >
>>> plt.show()

生成的图像如下: Output of GeoPandas nybb.plot() 这种快速绘图有助于检查数据是否正确加载以及几何图形是否合理。但通常情况下,您需要更进一步,使用颜色来显示每个区域的值。根据列中的值对区域进行着色的地图称为等值线图,您可以通过将要可视化的列传递给 column 参数来创建它:

>>> nybb.plot(column="Shape_Area", legend=True)
<Axes: >
>>> plt.show()

这将产生以下输出: Output of GeoPandas nybb.plot(column="Shape_Area", legend=True) 在这里,你可以看到每个行政区的相对大小,并用颜色表示。

创建交互式地图

您还可以使用 .explore() 方法创建交互式、可缩放的地图。.explore() 方法底层使用Folium生成 HTML 地图。REPL 无法直接显示该地图,因此 .show_in_browser() 方法会在您的默认浏览器中打开地图:

>>> nybb.explore().show_in_browser()

这将在您的浏览器中打开一个交互式地图。您应该会看到类似这样的内容: Output of GeoPandas .explore() Method on nybb Dataset 与绘制简单静态形状的.plot()不同,.explore()将您的数据叠加到参考图纸上——参考图纸包含地名、海岸线和缩放控件。这使得在上下文中解读数据以及交互式探索地理环境变得更加容易。

.plot() 类似,您可以通过传递列名将其转换为等值线图。这里,您可以根据每个行政区的土地面积为其着色:

>>> nybb.explore(column="Shape_Area", legend=False).show_in_browser()

输出结果为: GeoPandas .explore() Method on nybb Dataset 现在,每个行政区的相对大小用颜色表示,让您了解交互式地图如何揭示数据中的模式。为了使这些地图准确且在不同数据集之间具有可比性,接下来您将了解坐标参考系统和投影的工作原理。

使用坐标参考系统 (CRS) 和投影

每个 GeoDataFrame 都包含一个坐标参考系统 (CRS,Coordinate Reference Systems ) ,用于解释其坐标如何映射到现实世界的位置。理解 CRS 对于解释几何值、比较数据集和生成精确地图至关重要。在进一步探索之前,请仔细查看 nybb 中的 geometry 列:

>>> nybb[["BoroName", "geometry"]]
     BoroName       geometry
0    Staten Island  MULTIPOLYGON (((970217.022 145643.332, 970227....
1    Queens         MULTIPOLYGON (((1029606.077 156073.814, 102957...
2    Brooklyn       MULTIPOLYGON (((1021176.479 151374.797, 102100...
3    Manhattan      MULTIPOLYGON (((981219.056 188655.316, 980940....
4    Bronx          MULTIPOLYGON (((1012821.806 229228.265, 101278...

请注意,geometry 列只是数字——这些值如何与地球上的地点对应并不明显。

每个 GeoDataFrame 的 CRS 位于其 .crs 属性中。以下是 nybb 数据集的 CRS:

>>> nybb.crs
<Projected CRS: EPSG:2263>
Name: NAD83 / New York Long Island (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - New York - counties of Bronx; Kings; Nassau ...
- bounds: (-74.26, 40.47, -71.8, 41.3)
Coordinate Operation:
- name: SPCS83 New York Long Island zone (US survey foot)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

每个坐标参考系 (CRS) 都有一个名称和一个代码。在本例中,CRS 的名称为 NAD83 / New York Long Island (ftUS),代码为 EPSG:2263。x 和 y 坐标表示以美国测量英尺为单位的距离,原点位于纽约市郊外很远的地方。这种位置关系解释了为什么坐标轴的值很大,因为它们表示每个点距离原点的东向和北向距离(以英尺为单位)。

EPSG:2263 是一个投影坐标参考系。这意味着 geometry 列中的坐标已经过数学投影转换,使得 x 轴和 y 轴表示线性距离。下一节将介绍地理坐标参考系,其中坐标轴并非表示距离,而是表示纬度和经度。

探索地理 CRS

查看世界地图时常用的坐标参考系统 (CRS) 是 WGS 84 (EPSG:4326)。以下是如何使用此 CRS 加载世界地图的方法:

>>> world = gpd.read_file(get_path("naturalearth land"))
>>> world.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

WGS 84 是一种地理坐标参考系:x 轴和 y 轴分别代表经度和纬度(单位为度)。它们并不代表英尺或米等实际距离。由于度数并非均匀距离,因此直接使用地理坐标参考系绘制图形会扭曲形状和面积,尤其是在两极地区。

以下是绘制世界地图的方法:

>>> world.plot()
<Axes: >
>>> plt.show()

结果如下: A World Map Using the WGS84 Geographic Projection 由于该 CRS 不是投影的,因此极地附近的区域看起来被拉伸了——格陵兰岛和南极洲看起来比实际大得多。

改变 CRS

您可以更改 GeoDataFrame 的坐标参考系统 (CRS)。这很重要,因为 CRS 的选择会影响距离、面积和形状的显示方式。

例如,您刚才看到,WGS 84 投影使格陵兰岛和南极洲看起来比实际大得多。莫尔维德投影(Mollweide projection)是一种等面积投影:每个区域都保持其相对于其他区域的真实面积比例。形状变得更加圆润,距离也会有所变化,但各大洲和国家的相对大小是准确的。

您可以使用方法 .to_crs() 将世界地图的 CRS 从 WGS 84 更改为 Mollweide:

>>> world_moll = world.to_crs("ESRI:54009")
>>> world_moll.crs
<Projected CRS: ESRI:54009>
Name: World_Mollweide
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Coordinate Operation:
- name: World_Mollweide
- method: Mollweide
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

与所有投影式 CRS 一样,地图的坐标轴以距离表示,而不是以纬度和经度表示。

你可以像这样使用新的CRS绘制世界地图:

>>> world_moll.plot()
<Axes: >
>>> plt.show()

输出结果如下: Output of GeoPandas .plot() using a Mollweide Projection 在这张地图中,格陵兰岛和南极洲看起来比在WGS 84坐标系中小得多。更重要的是,它们的大小与其他陆地相比是准确的。

比较 CRS

有时,对两种不同坐标参考系下的数据集进行并排比较会很有帮助。以下示例使用Matplotlib ,分别使用 WGS 84 和 Mollweide 坐标系绘制了同一张世界地图:

>>> fig, axes = plt.subplots(1, 2, figsize=(14, 7))
>>> world.plot(ax=axes[0])
<Axes: >
>>> axes[0].set_title("EPSG:4326 (WGS84 Lat/Lon)")
Text(0.5, 1.0, 'EPSG:4326 (WGS84 Lat/Lon)')
>>> world_moll.plot(ax=axes[1])
<Axes: >
>>> axes[1].set_title("ESRI:54009 (Mollweide Equal-Area)")
Text(0.5, 1.0, 'ESRI:54009 (Mollweide Equal-Area)')
>>> plt.tight_layout()
>>> plt.show()

这就是最终的对比图: Comparison of two Projections of a World Map 从图中可以明显看出,坐标参考系统 (CRS) 的选择会影响地图的渲染方式。由于 WGS 84 将经度和纬度视为 x 坐标和 y 坐标,因此极地附近的区域看起来被拉伸,比实际面积更大。相比之下,莫尔维德投影保留了各陆地相对大小,因此格陵兰岛和南极洲看起来要小得多。

理解地图投影和中心参考系统

正如您所见,GeoPandas 使用**坐标参考系统 (CRS)**将坐标映射到现实世界的位置。有些 CRS 是地理坐标系,使用以度为单位的纬度和经度。这些坐标系未经投影。另一些坐标系则经过投影,这意味着它们使用地图投影将地球平面化为以线性单位表示的 x 和 y 坐标。

每个投影都是 CRS 的一部分,但并非每个 CRS 都使用投影。

执行空间连接(join)

空间连接允许您根据两个 GeoDataFrame 的几何关系合并它们的信息。与通过共享 ID 或列名匹配行不同,空间连接根据要素在空间中的位置进行匹配——例如,哪些地铁站位于某个行政区内。它相当于数据库连接的地理空间版本,但关键在于位置而非字段值。

在下面的示例中,您将使用空间连接来确定纽约市的哪个行政区包含帝国大厦。首先创建一个包含帝国大厦经度和纬度的 GeoDataFrame。由于这些坐标以度为单位,因此您应该使用 EPSG:4326 作为 CRS:

>>> from shapely.geometry import Point
>>> empire_state = gpd.GeoDataFrame(
...    [{"name": "Empire State Building"}],
...    geometry=[Point(-73.9857, 40.7484)],
...    crs="EPSG:4326"
... )
>>> empire_state
     name                     geometry
0    Empire State Building    POINT (-73.9857 40.7484)

现在您拥有一个只包含一个点的 GeoDataFrame。

理解空间连接和 CRS

空间连接要求两个 GeoDataFrame 具有相同的 CRS(中心参考系)。您刚才看到 empire_state DataFrame 使用了 EPSG:4326。之前您看到 nybb 使用了 EPSG:2263。以下是如何将 empire_state 也更改为使用 EPSG:2263

>>> empire_state = empire_state.to_crs(nybb.crs)
>>> empire_state
     name                     geometry
0    Empire State Building    POINT (988212.237 211939.279)

请注意,当您更改 empire_state 的坐标参考系 (CRS) 时,geometry 列中的值从 (-73.985740.7484) 变为 (988212.237211939.279)。这是因为坐标只有在特定的坐标参考系 (CRS) 下才有意义。当您切换到不同的坐标参考系 (CRS) 时,同一个实际位置会用不同的坐标值来表示。

使用.sjoin()方法

GeoPandas 中的空间连接与pandas 中的连接方式类似。具体来说,GeoDataFrame 有一个 .sjoin() 实例方法。此方法需要您指定要连接的 GeoDataFrame,以及要使用的连接类型和谓词

在下面的示例中,您使用左连接将 empire_state DataFrame 连接到 nybb DataFrame。使用左连接,您至少可以得到一行结果,该行将来自原始 empire_state DataFrame。

在运行连接操作之前,了解参数 predicate 的工作原理非常重要。参数 predicate 告诉 .sjoin() 如何判断两个几何图形在空间上是否匹配。这里,您使用 "within" 是因为您想测试一个点是否位于某个行政区内:

>>> empire_state.sjoin(nybb, how="left", predicate="within")
   name                   geometry               index_right  BoroName  ...
0  Empire State Building  POINT (988212.237 ...  3            Manhattan ...

由于纽约市各行政区互不重叠,所以只返回一行数据。帝国大厦位于曼哈顿,连接操作也显示了这一点。由于您连接了两个数据框,因此您同时拥有来自 empire_statenybb 的列。

在这个例子中,你了解了"within"谓词,它可以用来测试一个形状是否包含在另一个形状之内。其他谓词可以用来测试两个形状是否接触、重叠、相交等等。

避免局限性和陷阱

尽管 GeoPandas 功能强大且对初学者友好,但还是有一些常见的陷阱可能会让人犯错——尤其是在开始处理真实世界数据时。

记住这些要点可以节省您的时间和避免困惑:

  • **依赖项:**安装 GeoPandas 可能比较棘手,因为它依赖于已编译的库(GEOSPROJGDAL)。如果您在使用 pip 安装时遇到问题,请参阅官方安装页面以获取其他安装方案的信息。
  • **坐标参考系 (CRS) 注意:**务必检查并设置 CRS。在进行空间连接或距离计算之前忘记重新投影可能会导致错误的结果。
  • 地图绘制限制: .plot().explore() 非常适合快速绘制地图,但自定义功能有限。如需多图层、弹出窗口和底图等高级交互功能,建议直接使用 Folium。

这些都不是什么大问题,这些是值得及早识别的模式,这样你的工作流程才能保持流畅和可预测。

结论

GeoPandas 将熟悉的 pandas 工作流程扩展到了地理空间领域。在此过程中,您学习了如何加载常见的地理空间格式、创建静态和交互式地图、使用和更改数据集的 CRS,以及使用空间连接来回答基于位置的问题。

您现在已掌握处理真实世界空间数据并创建能够传达有意义见解的地图所需的核心技能。您可以通过尝试使用自己的数据集、阅读GeoPandas 文档或参考 Real Python 的Folium 教程来创建交互式 Web 地图,从而进一步提升这些技能。

realpython.com/geopandas/

mp.weixin.qq.com/s/-1MHFlFpf…