GeoPandas快速入门
官方文档:https://geopandas.org/en/stable/getting_started/introduction.html
GeoPandas 是一个开源项目,旨在为 pandas 对象添加对地理数据的支持。它当前实现 GeoSeries
和 GeoDataFrame
类型,它们分别是 pandas.Series
和 pandas.DataFrame
的子类。 GeoPandas 对象可以作用于 shapely
几何对象并执行几何操作。
一、概念
GeoPandas 通过添加对地理空间数据的支持来扩展流行的数据科学库 pandas。如果您不熟悉 pandas
,我们建议您在继续之前快速查看其入门文档。
GeoPandas 中的核心数据结构是 geopandas.GeoDataFrame
,它是 pandas.DataFrame
的子类,可以存储几何列并执行空间操作。geopandas.GeoSeries
是 pandas.Series
的子类,用于处理几何。
GeoDataFrame
是 pandas.Series
传统数据(数值、布尔值、文本等)和 geopandas.GeoSeries
与几何图形(点、多边形等)的组合。可以根据需要设置任意多个包含几何图形的列。
每个 GeoSeries
可以包含任何几何类型(甚至可以将它们混合在一个数组中)并具有 GeoSeries.crs
属性,该属性存储有关投影的信息(CRS 代表坐标参考系统) 。因此, GeoDataFrame
中的每个 GeoSeries
可以位于不同的投影中,例如,允许拥有同一几何图形的多个版本(不同的投影)。
GeoDataFrame
中只有一个 GeoSeries
被视为活动几何图形,这意味着应用于 GeoDataFrame
的所有几何操作都在此活动列上进行。
二、读取和写入文件
读取文件
假设您有一个包含数据和几何图形的文件(例如 GeoPackage、GeoJSON、Shapefile),您可以使用 geopandas.read_file()
读取它,它会自动检测文件类型并创建 GeoDataFrame
。本教程使用 "nybb"
数据集,这是纽约行政区的地图,可通过 geodatasets
包获取。因此,我们使用 geodatasets.get_path()
下载数据集并检索本地副本的路径。
conda install geodatasets //安装geodatasets
import geopandas
from geodatasets import get_path
path_to_data = get_path("nybb")
gdf = geopandas.read_file(path_to_data)
gdf
写入文件
要将 GeoDataFrame
写回文件,请使用 GeoDataFrame.to_file()
。默认文件格式是 Shapefile,但您可以使用 driver
关键字指定您自己的文件格式。
gdf.to_file("my_file.geojson", driver="GeoJSON")
简单的访问器和方法
现在我们有了 GeoDataFrame
并且可以开始处理它的几何形状。
由于 New York Boroughs 数据集中只有一个几何列,因此该列自动成为活动几何,并且 GeoDataFrame
上使用的空间方法将应用于 "geometry"
列。
测量面积
要测量每个多边形(在本例中为 MultiPolygon)的面积,请访问 GeoDataFrame.area
属性,该属性返回 pandas.Series
。请注意, GeoDataFrame.area
只是应用于活动几何列的 GeoSeries.area
。
但首先,为了使结果更易于阅读,请将行政区的名称设置为索引:
gdf = gdf.set_index("BoroName")
gdf["area"] = gdf.area
gdf["area"]
获取多边形边界和质心
要获取每个多边形的边界 (LineString),使用GeoDataFrame.boundary
:
gdf["boundary"] = gdf.boundary
gdf["boundary"]
由于我们已将边界保存为新列,因此现在GeoDataFrame
中有两个几何列。
我们还可以创建新的几何图形,例如,它可以是原始几何图形的缓冲版本(即, GeoDataFrame.buffer(10)
)或其质心:
gdf["centroid"] = gdf.centroid
gdf["centroid"]
测量距离
测量每个质心距离第一个质心位置有多远。
first_point = gdf["centroid"].iloc[0]
gdf["distance"] = gdf["centroid"].distance(first_point)
gdf["distance"]
请注意, geopandas.GeoDataFrame
是 pandas.DataFrame
的子类,因此我们可以在地理空间数据集上使用所有pandas功能-我们甚至可以同时使用属性和几何信息执行数据操作。
例如,要计算上面测量的距离的平均值,调用“distance”列并调用其mean()
方法:
gdf["distance"].mean()
制作地图
GeoPandas 还可以绘制地图,因此我们可以查看几何对象在空间中的外观。要绘制活动几何对象,请调用 GeoDataFrame.plot()
。要按另一列进行颜色编码,请将该列作为第一个参数传递。在下面的示例中,我们绘制活动几何列,并按"area
"列进行颜色编码。我们还想显示一个图例(legend=True
)。
gdf.plot("area", legend=True)
还可以使用 GeoDataFrame.explore()
交互地探索数据,其行为方式与 plot()
相同,但返回交互式地图。
gdf.explore("area", legend=False)
遇到问题:报错:CRSError: Invalid projection: EPSG:4326: (Internal Proj Error: proj_create: no database context specified)
解决方法[^1]
将活动的几何( GeoDataFrame.set_geometry
)切换到质心,我们可以使用点绘制相同的数据。
gdf = gdf.set_geometry("centroid")
gdf.plot("area", legend=True)
我们也可以将两者叠加在一起。我们只需要用一个图作为另一个图的轴。
ax = gdf["geometry"].plot()
gdf["centroid"].plot(ax=ax, color="black")
现在,我们将活动几何体设置回原始的 GeoSeries
。
gdf = gdf.set_geometry("geometry")
几何创建
凸包
如果我们对多边形的凸包感兴趣,我们可以访问 GeoDataFrame.convex_hull
。
gdf["convex_hull"] = gdf.convex_hull
# saving the first plot as an axis and setting alpha (transparency) to 0.5
ax = gdf["convex_hull"].plot(alpha=0.5)
# passing the first plot and setting linewidth to 0.5
gdf["boundary"].plot(ax=ax, color="white", linewidth=0.5)
缓冲
在其他情况下,我们可能需要使用 GeoDataFrame.buffer()
对几何图形进行缓冲。几何方法会自动应用于活动几何图形,但我们也可以直接将它们应用于任何 GeoSeries
。让我们对行政区和它们的质心进行缓冲,并将它们叠加在一起绘制出来。
# buffering the active geometry by 10 000 feet (geometry is already in feet)
gdf["buffered"] = gdf.buffer(10000)
# buffering the centroid geometry by 10 000 feet (geometry is already in feet)
gdf["buffered_centroid"] = gdf["centroid"].buffer(10000)
# saving the first plot as an axis and setting alpha (transparency) to 0.5
ax = gdf["buffered"].plot(alpha=0.5)
# passing the first plot as an axis to the second
gdf["buffered_centroid"].plot(ax=ax, color="red", alpha=0.5)
# passing the first plot and setting linewidth to 0.5
gdf["boundary"].plot(ax=ax, color="white", linewidth=0.5)
几何关系
我们还可以询问不同几何形状的空间关系。使用上面的几何图形,我们可以检查哪些缓冲行政区与布鲁克林的原始几何图形相交,即距布鲁克林 10000 英尺以内。
首先,我们得到布鲁克林的多边形。
brooklyn = gdf.loc["Brooklyn", "geometry"]
brooklyn
这个多边形是一个 Shapely 几何对象,就像 GeoPandas 中使用的其他几何对象一样。
type(brooklyn)
# shapely.geometry.multipolygon.MultiPolygon
然后我们可以检查 gdf["buffered"]
中的哪个几何图形与它相交。
gdf["buffered"].intersects(brooklyn)
只有布朗克斯(在北部)距离布鲁克林超过 10,000 英尺。其他所有地区都更接近并与我们的多边形相交。
或者,我们可以检查哪些缓冲质心完全位于原始行政区多边形内。在这种情况下,两个 GeoSeries
都是对齐的,并且对每一行执行检查。
gdf["within"] = gdf["buffered_centroid"].within(gdf)
gdf["within"]
我们可以将结果绘制在地图上.
gdf = gdf.set_geometry("buffered_centroid")
# using categorical plot and setting the position of the legend
ax = gdf.plot(
"within", legend=True, categorical=True, legend_kwds={"loc": "upper left"}
)
# passing the first plot and setting linewidth to 0.5
gdf["boundary"].plot(ax=ax, color="black", linewidth=0.5)
投影
每个 GeoSeries
都有其坐标参考系统 (CRS),可通过 GeoSeries.crs
访问。 CRS 告诉 GeoPandas 几何图形的坐标位于地球表面的位置。在某些情况下,CRS 是地理的,这意味着坐标采用纬度和经度。在这些情况下,其 CRS 为 WGS84,权限代码为 EPSG:4326
。让我们看看纽约行政区 GeoDataFrame
的投影。
gdf.crs
几何图形位于 EPSG:2263
中,坐标以英尺为单位。我们可以轻松地将 GeoSeries
重新投影到另一个 CRS,例如使用 GeoSeries.to_crs()
的 EPSG:4326
。
gdf = gdf.set_geometry("geometry")
boroughs_4326 = gdf.to_crs("EPSG:4326")
boroughs_4326.plot()
boroughs_4326.crs
请注意绘图坐标轴上的坐标差异。之前我们是以 120,000 到 280,000(英尺)为单位,现在是以 40.5 到 40.9(度)为单位。在这种情况下,boroughs_4326
具有 WGS84 中的“geometry
”列,但所有其他列(包括质心等)仍保持在原始 CRS 中。
更多
借助 GeoPandas,我们可以做的事情比迄今为止介绍的要多得多,从聚合到空间连接,再到地理编码等等。
用户手册:https://geopandas.org/en/stable/docs/user_guide.html