【GIS】GeoPandas：Python矢量数据处理

GeoPandas 是 Python 中处理矢量地理数据的强大工具。它结合了 Pandas 的数据处理能力和 Shapely 几何对象的处理功能，使得空间数据的处理和分析变得简单而高效。本文将详细介绍如何使用 GeoPandas 创建数据、操作几何对象、进行空间查询以及可视化展示，帮助您快速上手并掌握地理数据处理。

数据结构

以下是对 Point、LineString、Polygon、DataFrame 和 GeoDataFrame 的简单介绍，包括它们的内置方法和常见用法。

Point

Point 表示二维空间中的一个坐标点。

from shapely.geometry import Point# 创建一个点对象
point = Point(1, 2)

常用方法:

• x: 获取点的 x 坐标。
• y: 获取点的 y 坐标。
• distance(other): 计算当前点到另一个点的距离。
• buffer(radius): 创建一个以当前点为中心的缓冲区。

LineString

LineString 表示由一系列点连接成的线段。

from shapely.geometry import LineString# 创建一个线对象
line = LineString([(0, 0), (1, 2), (2, 1)])

常用方法:

• length: 获取线的长度。
• coords: 获取线的坐标点列表。
• intersects(other): 判断当前线是否与其他几何对象相交。
• buffer(radius): 创建线的缓冲区。

Polygon

Polygon 表示由多个点围成的多边形。

from shapely.geometry import Polygon# 创建一个多边形对象
polygon = Polygon([(0, 0), (4, 0), (4, 3), (0, 3)])

常用方法:

• area: 获取多边形的面积。
• perimeter: 获取多边形的周长。
• contains(other): 判断多边形是否包含另一个几何对象。
• intersection(other): 计算与另一个几何对象的交集。

DataFrame

DataFrame 是 Pandas 中的核心数据结构，用于存储表格数据。

import pandas as pd# 创建一个简单的 DataFrame
data = {'Name': ['A', 'B'], 'Value': [10, 20]}
df = pd.DataFrame(data)

常用方法:

• head(n): 返回前 n 行数据。
• describe(): 生成描述性统计信息。
• groupby(by): 按指定列进行分组。
• filter(): 根据条件筛选数据。

GeoDataFrame

GeoDataFrame 是 GeoPandas 中的扩展数据结构，专门用于存储地理数据。

import geopandas as gpd# 创建 GeoDataFrame
gdf = gpd.GeoDataFrame({'geometry': [point, line, polygon]})

常用方法:

• set_crs(crs): 设置坐标参考系统（CRS）。
• to_crs(crs): 转换坐标参考系统。
• sjoin(other, op): 进行空间连接。
• plot(): 可视化 GeoDataFrame。

地理数据操作

创建新数据

您可以通过手动创建几何对象和相关属性来新建一个 GeoDataFrame。以下是一个简单的示例：

import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon# 创建几何对象
point = Point(1, 2)  # 一个点
line = LineString([(0, 0), (1, 2), (2, 1)])  # 一条线
polygon = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])  # 一个多边形# 创建 GeoDataFrame
gdf = gpd.GeoDataFrame({'geometry': [point, line, polygon]})# 打印 GeoDataFrame
print(gdf)

从地理数据导入

您可以从文件中读取地理数据，例如 Shapefile 或 GeoJSON 等格式。

# 导入 GeoDataFrame
gdf = gpd.read_file("path/to/your/file.shp")
print(gdf.head())

从表格导入

如果您的数据存储在 Excel 或 CSV 文件中，可以轻松导入并转换为 GeoDataFrame。

import pandas as pd# 从 CSV 文件读取数据
df = pd.read_excel("path/to/your/file.csv")# 创建 GeoDataFrame，指定经纬度列
point = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lng, df.lat))
point.crs = 'EPSG:4326'  # 设置坐标系

导出数据

处理完成后，您可以将 GeoDataFrame 导出为多种地理数据格式。

# 导出 GeoDataFrame
gdf.to_file("path/to/save/file.geojson", driver='GeoJSON')

数据查看

查看数据

使用描述性统计和可视化方法对数据进行探索。

# 查看数据的前几行
print(gdf.head())# 查看数据的列名
print(gdf.columns)# 查看数据的几何类型
print(gdf.geom_type)

绘制基本图形

您可以使用 GeoPandas 内置的绘图功能来可视化数据。

import matplotlib.pyplot as plt# 绘制 GeoDataFrame
gdf.plot()
plt.title("Basic GeoDataFrame Plot")
plt.show()

自定义可视化

您可以根据需要自定义绘图样式，以增强可视化效果。

# 自定义绘图样式
gdf.plot(color='lightblue', edgecolor='black')
plt.title("Custom GeoDataFrame Plot")
plt.show()

空间分析

缓冲区

您可以生成几何对象的缓冲区，以创建一定范围的区域。

# 生成缓冲区
gdf['buffer'] = gdf.geometry.buffer(0.5)# 可视化缓冲区
gdf.set_geometry('buffer').plot()
plt.title("Buffer Zones")
plt.show()

几何操作

创建两个示例 GeoDataFrame。

# 创建示例多边形
poly1 = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])  # 多边形 1
poly2 = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])  # 多边形 2# 创建 GeoDataFrame
gdf1 = gpd.GeoDataFrame({'geometry': [poly1]})
gdf2 = gpd.GeoDataFrame({'geometry': [poly2]})

使用 gpd.overlay 函数进行不同的几何操作，只需更改 how 参数：

• 交集（Intersection）:

  result = gpd.overlay(gdf1, gdf2, how='intersection')
  

• 并集（Union）:

  result = gpd.overlay(gdf1, gdf2, how='union')
  

• 异或（Symmetrical Difference）:

  result = gpd.overlay(gdf1, gdf2, how='symmetric_difference')
  

• 剪去（Difference）:

  result = gpd.overlay(gdf1, gdf2, how='difference')
  

可视化结果的方法相同，可以使用以下代码展示结果：

result.plot()
plt.title("Result Title")  # 根据具体操作修改标题
plt.show()

融合

您可以将多个几何对象融合为一个复合几何对象，这在处理相邻区域时特别有用。

# 创建多个几何对象
polygons = [Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),Polygon([(1, 0), (2, 0), (2, 1), (1, 1)])]# 创建 GeoDataFrame
gdf_polygons = gpd.GeoDataFrame({'geometry': polygons})# 融合几何对象
merged = gdf_polygons.unary_union# 可视化融合结果
gpd.GeoSeries([merged]).plot()
plt.title("Merged Geometry")
plt.show()

空间查询

除了通过 gpd.overlay 进行的几何操作，GeoPandas 还支持多种空间查询操作，包括：

• 相交查询（Intersects）:
  检查哪些几何对象相交。

  intersects = gdf1[gdf1.geometry.intersects(gdf2.geometry.iloc[0])]
  

• 包含查询（Contains）:
  检查哪些几何对象包含其他几何对象。

  contains = gdf1[gdf1.geometry.contains(gdf2.geometry.iloc[0])]
  

• 相离查询（Does Not Intersect）:
  检查哪些几何对象不相交。

  disjoint = gdf1[gdf1.geometry.disjoint(gdf2.geometry.iloc[0])]
  

• 重叠查询（Overlaps）:
  检查哪些几何对象部分重叠。

  overlaps = gdf1[gdf1.geometry.overlaps(gdf2.geometry.iloc[0])]
  

• 触碰查询（Touches）:
  检查哪些几何对象相接触但不重叠。

  touches = gdf1[gdf1.geometry.touches(gdf2.geometry.iloc[0])]
  

对于每种查询，您可以使用类似的可视化方法：

intersects.plot(color='blue')
plt.title("Intersecting Geometries")
plt.show()

空间连接

在使用 sjoin 进行空间连接后，返回的 GeoDataFrame 会包含两个 GeoDataFrame 的所有列，同时会增加一个用于表示空间关系的列（例如，index_right，指向右侧 GeoDataFrame 的索引）。以下是 sjoin 的基本使用示例和结果分析。

import geopandas as gpd
from shapely.geometry import Point, Polygon# 创建第一个 GeoDataFrame
gdf1 = gpd.GeoDataFrame({'name': ['A', 'B', 'C'],'geometry': [Point(1, 1), Point(2, 2), Point(3, 3)]
})# 创建第二个 GeoDataFrame
gdf2 = gpd.GeoDataFrame({'name': ['Zone 1', 'Zone 2'],'geometry': [Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]), Polygon([(2, 2), (4, 2), (4, 4), (2, 4)])]
})# 进行空间连接
joined = gpd.sjoin(gdf1, gdf2, how='inner', op='intersects')# 打印结果
print(joined)

在上述代码中，sjoin 会返回一个新的 GeoDataFrame，结果可能类似于：

  name_left     geometry_left name_right                  geometry_right
0        A      POINT (1 1)     Zone 1  POLYGON ((0 0, 2 0, 2 2, 0 2))
1        B      POINT (2 2)     Zone 1  POLYGON ((0 0, 2 0, 2 2, 0 2))
2        C      POINT (3 3)     Zone 2  POLYGON ((2 2, 4 2, 4 4, 2 4))

• 列名: name_left 和 name_right 分别表示左侧和右侧 GeoDataFrame 的列。
• 几何对象: geometry_left 和 geometry_right 分别表示两个 GeoDataFrame 中的几何对象。
• 索引: index_right 列指向与左侧对象相交的右侧 GeoDataFrame 的索引。

坐标系操作

GeoPandas 提供了方便的坐标系操作功能，帮助您轻松管理和转换不同的地理数据坐标系。

定义与转换坐标系

您可以为 GeoDataFrame 定义坐标参考系统 (CRS)，并进行转换。

# 定义 CRS
gdf.crs = "EPSG:4326"  # WGS84# 转换为另一 CRS
gdf_transformed = gdf.to_crs("EPSG:3857")  # Web Mercator
print(gdf_transformed.head())

处理坐标系问题

确保在进行空间操作前，所有 GeoDataFrame 使用相同的 CRS，以避免潜在问题。

# 检查 CRS
print(gdf.crs)
print(gdf_transformed.crs)

常用坐标系

在 GeoPandas 中，常用坐标系通过 EPSG（欧洲规范化地理信息系统）代码来表示。以下是一些常用坐标系及其在 GeoPandas 中的表示方式：

坐标系名称	EPSG 代码	描述
WGS 84	4326	全球通用坐标系，使用经纬度表示。
Web Mercator	3857	常用于网络地图（如 Google Maps），单位为米。
CGCS 2000	4490	中国地理坐标系统。
UTM Zone 33N	32633	通用横轴墨卡托投影，适用于小范围区域（例如欧洲）。
国家测绘局坐标系	4547	用于中国地理数据，特别是在一些国家级项目中。
BD-09	N/A	百度地图使用的坐标系，无法直接用 EPSG 表示。