Python与GIS：使用GeoPandas和Shapely进行地理空间数据处理和分析
大家好！今天我们来深入探讨如何使用Python生态系统中的两个强大工具——GeoPandas和Shapely，进行地理空间数据的处理和分析。GIS（Geographic Information System，地理信息系统）在现代数据科学中扮演着越来越重要的角色，而Python以其易用性和丰富的库支持，成为了GIS分析的首选语言。
一、GIS数据基础与概念
在开始编码之前，我们需要了解一些基本的GIS概念：
-
几何对象（Geometry Objects）： 这是GIS数据的基本构建块，用于表示地球上的各种地理实体。常见的几何对象包括：
- 点（Point）： 表示一个单独的位置，例如城市中心或测量站。
- 线（LineString）： 表示一系列有序的点，连接起来形成一条线，例如道路或河流。
- 面（Polygon）： 表示一个封闭的区域，由一系列线段构成，例如国家边界或湖泊。
- 多点（MultiPoint）： 表示多个点的集合。
- 多线（MultiLineString）： 表示多条线的集合。
- 多面（MultiPolygon）： 表示多个面的集合。
-
坐标参考系统（Coordinate Reference System, CRS）： 用于定义地球表面上的位置。常见的CRS包括：
- 地理坐标系（Geographic CRS）： 使用经度和纬度来表示位置，例如WGS 84 (EPSG:4326)。
- 投影坐标系（Projected CRS）： 将地球表面投影到平面上，使用米或英尺作为单位，例如UTM (Universal Transverse Mercator)。
-
属性数据（Attribute Data）： 与几何对象关联的数据，用于描述地理实体的特征。例如，一个城市（Polygon）的属性数据可能包括人口、面积和行政级别。
二、Shapely：几何对象的操作利器
Shapely是一个Python库，专门用于操作和分析几何对象。它提供了创建、操作和分析点、线和面的功能。
1. 安装Shapely
pip install shapely2. 创建几何对象
from shapely.geometry import Point, LineString, Polygon
# 创建一个点
point = Point(2.5, 4.5)
print(point) # 输出: POINT (2.5 4.5)
print(point.x, point.y) # 输出 2.5 4.5
# 创建一条线
line = LineString([(0, 0), (1, 1), (2, 0), (3, 1), (4, 0)])
print(line) # 输出: LINESTRING (0 0, 1 1, 2 0, 3 1, 4 0)
# 创建一个面
polygon = Polygon([(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)])
print(polygon) # 输出: POLYGON ((0 0, 0 2, 2 2, 2 0, 0 0))3. 几何对象的属性和方法
Shapely提供了丰富的属性和方法来获取几何对象的各种信息和进行各种操作。
| 属性/方法 | 描述 |
|---|---|
area |
返回面的面积。 |
bounds |
返回几何对象的边界框 (minx, miny, maxx, maxy)。 |
centroid |
返回几何对象的质心。 |
distance(other) |
返回几何对象与另一个几何对象之间的最小距离。 |
length |
返回线的长度。 |
contains(other) |
如果几何对象包含另一个几何对象，则返回 True。 |
intersects(other) |
如果几何对象与另一个几何对象相交，则返回 True。 |
union(other) |
返回几何对象与另一个几何对象的并集。 |
intersection(other) |
返回几何对象与另一个几何对象的交集。 |
difference(other) |
返回几何对象与另一个几何对象的差集。 |
buffer(distance) |
返回几何对象的缓冲区，创建一个指定距离的缓冲区多边形。 |
simplify(tolerance) |
返回几何对象的简化版本，减少顶点数量，同时保持形状的近似。 tolerance 控制简化程度。 |
wkt |
返回几何对象的 WKT (Well-Known Text) 表示。 |
wkb |
返回几何对象的 WKB (Well-Known Binary) 表示。 |
# 示例
polygon = Polygon([(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)])
print(f"Area: {polygon.area}")
print(f"Bounds: {polygon.bounds}")
print(f"Centroid: {polygon.centroid}")
point = Point(1, 1)
print(f"Polygon contains point: {polygon.contains(point)}")
# 创建另一个多边形
polygon2 = Polygon([(1, 1), (1, 3), (3, 3), (3, 1), (1, 1)])
# 计算两个多边形的交集
intersection = polygon.intersection(polygon2)
print(f"Intersection: {intersection}")三、GeoPandas：地理空间数据分析的瑞士军刀
GeoPandas是建立在Pandas之上的一个库，它扩展了Pandas的数据结构，使其能够处理地理空间数据。GeoPandas的核心数据结构是GeoDataFrame，它类似于Pandas的DataFrame，但多了一列用于存储几何对象。
1. 安装GeoPandas
pip install geopandas2. 从文件读取地理空间数据
GeoPandas可以从多种地理空间数据格式读取数据，包括Shapefile、GeoJSON、PostGIS等。
import geopandas
# 从Shapefile读取数据
gdf = geopandas.read_file("path/to/your/shapefile.shp")
print(gdf.head())
print(gdf.crs) # 查看坐标参考系请将 "path/to/your/shapefile.shp" 替换为你实际的Shapefile路径。
3. GeoDataFrame的基本操作
GeoDataFrame继承了Pandas DataFrame的所有功能，并添加了地理空间数据特有的功能。
- 访问几何对象列：
# 访问几何对象列
geometry = gdf['geometry']
print(geometry.head())- 坐标参考系转换：
# 转换坐标参考系
gdf_wgs84 = gdf.to_crs("EPSG:4326") # 转换为WGS 84
print(gdf_wgs84.crs)- 空间连接（Spatial Join）： 根据空间关系将两个
GeoDataFrame连接起来。
# 假设我们有两个GeoDataFrame: gdf1 和 gdf2
# 进行空间连接，找到 gdf1 中每个要素所在的 gdf2 的要素
joined = geopandas.sjoin(gdf1, gdf2, how="left", predicate="intersects") # 或者 "contains", "within" 等
print(joined.head())sjoin函数的how参数指定连接方式（"left", "right", "inner"），predicate参数指定空间关系（"intersects", "contains", "within", "touches", "crosses", "overlaps"）。
- 空间分析： 使用
GeoPandas和Shapely进行各种空间分析操作。
# 计算每个几何对象的面积
gdf['area'] = gdf.geometry.area
# 创建缓冲区
gdf['buffer'] = gdf.geometry.buffer(1000) # 单位取决于坐标参考系
# 计算两个GeoDataFrame之间的空间关系
# 例如，找到所有与另一个GeoDataFrame相交的要素
intersecting = gdf[gdf.intersects(another_gdf.unary_union)]- 数据可视化：
GeoPandas可以方便地将地理空间数据可视化。
# 使用matplotlib进行可视化
import matplotlib.pyplot as plt
gdf.plot()
plt.show()
# 添加颜色和图例
gdf.plot(column='population', legend=True)
plt.show()4. 实战案例：分析城市人口分布
假设我们有两个Shapefile：一个是城市边界（city_boundary.shp），另一个是人口普查区（census_tracts.shp），每个普查区都有人口数据。我们想要计算每个城市的人口总数。
import geopandas
# 读取数据
city_boundary = geopandas.read_file("city_boundary.shp")
census_tracts = geopandas.read_file("census_tracts.shp")
# 确保两个GeoDataFrame使用相同的坐标参考系
census_tracts = census_tracts.to_crs(city_boundary.crs)
# 空间连接：将普查区连接到城市边界
joined = geopandas.sjoin(census_tracts, city_boundary, how="inner", predicate="intersects")
# 按城市名称分组，并计算人口总数
city_population = joined.groupby('city_name')['population'].sum()
print(city_population)
# 将结果添加到城市边界GeoDataFrame中
city_boundary['population'] = city_boundary['city_name'].map(city_population)
# 可视化
city_boundary.plot(column='population', legend=True)
plt.show()四、高级应用：地理空间数据处理的更多可能性
除了上述基本操作之外，GeoPandas和Shapely还可以用于更高级的地理空间数据处理和分析，例如：
- 网络分析： 使用
NetworkX库结合GeoPandas进行道路网络分析，例如最短路径计算、服务区分析等。 - 栅格数据处理： 使用
Rasterio库读取和处理栅格数据（例如卫星图像、DEM数据），并与GeoPandas矢量数据进行集成分析。 - 地理编码和反地理编码： 使用
Geopy库将地址转换为地理坐标（地理编码），或将地理坐标转换为地址（反地理编码）。 - 空间统计分析： 使用
PySAL库进行空间自相关分析、热点分析等。
五、一些常见问题的处理
使用GeoPandas和Shapely时，可能会遇到以下问题：
- 坐标参考系不匹配： 确保所有
GeoDataFrame使用相同的坐标参考系，否则空间分析结果可能不准确。使用to_crs()方法进行坐标参考系转换。 - 几何对象无效： 有时，几何对象可能由于数据错误或操作不当而变得无效。可以使用
Shapely的is_valid属性检查几何对象是否有效，并使用shapely.validation.make_valid()方法尝试修复无效几何对象。 - 空间连接性能问题： 对于大型数据集，空间连接可能非常耗时。可以尝试使用
Dask-GeoPandas库进行并行计算，提高空间连接的性能。 - 数据类型不一致： 在进行空间分析之前，确保所有数据类型都是正确的。例如，如果人口数据被错误地存储为字符串类型，需要将其转换为数值类型。
六、示例代码：计算两个GeoDataFrame的交集面积
import geopandas as gpd
from shapely.geometry import Polygon
# 创建两个示例GeoDataFrame
# 第一个GeoDataFrame
data1 = {'id': [1, 2],
'geometry': [Polygon([(0, 0), (0, 2), (2, 2), (2, 0)]),
Polygon([(1, 1), (1, 3), (3, 3), (3, 1)])]}
gdf1 = gpd.GeoDataFrame(data1, crs="EPSG:4326")
# 第二个GeoDataFrame
data2 = {'name': ['A', 'B'],
'geometry': [Polygon([(1, 0), (1, 2), (3, 2), (3, 0)]),
Polygon([(2, 1), (2, 3), (4, 3), (4, 1)])]}
gdf2 = gpd.GeoDataFrame(data2, crs="EPSG:4326")
# 确保两个GeoDataFrame使用相同的坐标参考系（这里已经相同了）
# 如果不相同，使用gdf2 = gdf2.to_crs(gdf1.crs)
# 定义一个函数来计算两个GeoDataFrame中每个要素的交集面积
def calculate_intersection_area(row):
return row['geometry'].intersection(gdf2.at[row['name'], 'geometry']).area
# 将gdf2的geometry列进行转换，name列作为索引
gdf2 = gdf2.set_index('name')
# 计算交集面积
# 创建一个包含所有组合的表
import pandas as pd
cross_join = pd.DataFrame([(i,j) for i in gdf1['id'] for j in gdf2.index], columns = ['id', 'name'])
cross_join = cross_join.merge(gdf1, on='id').merge(gdf2, left_on='name', right_index=True)
cross_join['intersection_area'] = cross_join.apply(calculate_intersection_area, axis=1)
print(cross_join)七、使用geopandas.clip()函数进行裁剪
import geopandas as gpd
from shapely.geometry import Polygon
# 创建一个示例GeoDataFrame (源数据)
data = {'id': [1, 2],
'geometry': [Polygon([(0, 0), (0, 4), (4, 4), (4, 0)]),
Polygon([(1, 1), (1, 3), (3, 3), (3, 1)])]}
gdf = gpd.GeoDataFrame(data, crs="EPSG:4326")
# 创建一个裁剪多边形 (裁剪区域)
clip_polygon = Polygon([(1, 1), (1, 3), (5, 3), (5, 1)])
clip_gdf = gpd.GeoDataFrame({'geometry': [clip_polygon]}, crs="EPSG:4326")
# 使用 geopandas.clip() 函数进行裁剪
clipped_gdf = gpd.clip(gdf, clip_polygon)
print(clipped_gdf)
# 可视化结果
import matplotlib.pyplot as plt
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))
gdf.plot(ax=ax1, color='blue', label='Original')
clip_gdf.plot(ax=ax1, color='red', alpha=0.5, label='Clip Polygon')
ax1.set_title('Original Data and Clip Polygon')
ax1.legend()
gdf.plot(ax=ax2, color='blue')
ax2.set_title('Original Data')
clipped_gdf.plot(ax=ax3, color='green')
ax3.set_title('Clipped Data')
plt.show()这个例子展示了如何使用 geopandas.clip() 函数将一个 GeoDataFrame 裁剪到指定的多边形区域。clip() 函数会返回一个新的 GeoDataFrame，其中包含原始 GeoDataFrame 中与裁剪多边形相交的部分。未相交的部分会被丢弃。
八、将结果保存到文件
# 将GeoDataFrame保存到Shapefile
gdf.to_file("output.shp", driver="ESRI Shapefile")
# 将GeoDataFrame保存到GeoJSON
gdf.to_file("output.geojson", driver="GeoJSON")确保安装了相应的驱动程序，例如，保存Shapefile需要fiona库。
九、灵活运用，掌握GIS分析的精髓
总而言之，GeoPandas和Shapely是Python中进行地理空间数据处理和分析的强大工具。通过灵活运用这些库，我们可以完成各种复杂的GIS任务，例如空间分析、数据可视化、网络分析等。希望今天的讲解能够帮助大家入门GeoPandas和Shapely，并在实际项目中应用它们。
十、空间数据处理的强大能力
GeoPandas和Shapely提供了处理空间数据的各种工具，可以进行创建、操作和分析几何对象。它们是Python中进行地理空间数据处理和分析的强大工具，可以完成各种复杂的GIS任务。掌握它们，能够使您灵活运用这些库，并在实际项目中应用它们。