地理空间数据库GeoPandas

概要

地理空间数据分析在许多领域中都有着广泛的应用，包括地理信息系统（GIS）、城市规划、环境科学、地质勘探等。Python GeoPandas是一个强大的库，旨在简化地理空间数据的处理和分析。本文将深入介绍GeoPandas库，包括其基本概念、功能特性、示例代码以及在实际应用中的应用场景。

一、什么是GeoPandas？

GeoPandas是一个开源的Python库，用于处理地理空间数据。它结合了两个其他流行的地理空间库，即Pandas和Shapely，提供了一个统一的数据结构来处理地理空间数据。GeoPandas的核心数据结构是GeoDataFrame，它是Pandas DataFrame的扩展，允许存储地理空间几何图形和属性数据。geopandas主要结合了pandas和shapely框架的能力。

  • shapely 有一个名为 geometry 的类,也是python 中一个非常重要的类库，其中包含不同的几何对象。
  • geopandas两个主要的数据结构，分别是GeoSeries和GeoDataFrame，分别对应了pandas中Series和DataFrame的子类。
  • geopandas有三个基本类的几何对象（实际上是形状对象）：点/点集合、 线/线集合、 多边形/多边形集合。

GeoPandas的主要功能包括：

  • 读取和写入多种地理空间数据格式，如Shapefile、GeoJSON、PostGIS、KML等。
  • 支持地理空间几何图形的创建、编辑和分析，包括点、线、面等。
  • 提供空间连接、缓冲区分析、几何运算等地理空间操作。
  • 具备数据可视化能力，可以绘制地理空间数据的地图。

二、安装GeoPandas

要开始使用GeoPandas，首先需要安装它。可以使用pip来安装GeoPandas及其依赖项：

pip install geopandas

安装完成后，可以在Python项目中引入GeoPandas并开始使用。

三、基本概念

3.1.GeoDataFrame

GeoPandas的核心数据结构是GeoDataFrame，它类似于Pandas的DataFrame，但包括一个额外的"geometry"列，用于存储地理空间几何图形。一个GeoDataFrame可以包含多个地理要素，每个要素具有一个几何图形和一组属性。以下是一个示例，展示了如何创建一个GeoDataFrame：

import geopandas as gpd
from shapely.geometry import Point

# 创建一个空的GeoDataFrame并指定几何列
gdf = gpd.GeoDataFrame(geometry=gpd.GeoSeries())

# 添加Point几何对象和属性数据
gdf['geometry'] = [Point(0, 0), Point(1, 1)]
gdf['name'] = ['Point A', 'Point B']

# 显式设置几何列,显式设置geometry列为GeoDataFrame的几何列，以确保GeoDataFrame正确识别几何数据。
gdf = gdf.set_geometry('geometry')

print(gdf)

上面的代码演示了如何使用GeoPandas库来创建和操作地理空间数据。具体来说，它展示了如何创建一个GeoDataFrame，向其中添加点几何对象和相关属性，并设置几何列。

3.2.几何图形

GeoPandas支持各种地理空间几何图形，包括点、线、面、多边形等。这些几何图形可以用来表示地理空间对象，如城市、河流、国家边界等。以下是一个示例，展示了如何创建不同类型的几何图形：

import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon

# 创建几何对象
"""
(1, 1)：这个点的坐标。通常，坐标是以 (x, y) 的形式表示，其中 x 是横坐标，y 是纵坐标。在地理空间数据中，x 和 y 分别对应于经度和纬度，或者在平面坐标系中对应于水平和垂直位置。
意义：表示一个具体的位置。例如，一个城市的具体位置、一个传感器的放置点等
"""
point = Point(1, 1)  # 创建一个点
"""
[(2, 0), (2, 4), (3, 4)]：这是一组有序的坐标点列表。每个点的坐标以 (x, y) 的形式表示。
意义：这些点按顺序连接形成一条线。表示线性地理对象，例如河流、道路、轨道等。
"""
line = LineString([(2, 0), (2, 4), (3, 4)])  # 创建一条线
"""
[(0, 0), (1, 0), (1, 1), (0, 1)]：这是一组有序的坐标点列表，定义多边形的顶点。第一个点和最后一个点通常是相同的，以闭合多边形。
意义：这些点按顺序连接形成一个闭合的区域。表示区域性地理对象，例如湖泊、国家边界、建筑轮廓等。
"""
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])  # 创建一个多边形

# 创建一个包含几何对象的GeoDataFrame
gdf = gpd.GeoDataFrame(geometry=[point, line, polygon])

# 添加属性数据
gdf['name'] = ['Point', 'Line', 'Polygon']

# 输出GeoDataFrame以查看其内容
print(gdf)

这些坐标数据是几何对象的基础，定义了对象在空间中的位置和形状：

  • 点 (Point)：单个坐标点，表示一个具体位置。
  • 线 (LineString)：由一系列坐标点按顺序连接形成的线，表示线性地理对象。
  • 多边形 (Polygon)：由一系列坐标点按顺序连接形成的闭合区域，表示面状地理对象

3.3.投影

地理空间数据通常需要投影到平面坐标系，以进行空间分析。GeoPandas支持各种地图投影，可以在需要时进行投影变换。以下是一个示例，展示了如何将GeoDataFrame投影到指定的投影坐标系：

import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon

# 创建几何对象
point = Point(1, 1)  # 创建一个点对象，坐标为(1, 1)
line = LineString([(2, 0), (2, 4), (3, 4)])  # 创建一个线对象，包含一系列坐标点
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])  # 创建一个多边形对象，包含一系列坐标点，形成一个闭合区域

# 创建一个包含几何对象的GeoDataFrame，并指定初始坐标系 (WGS84, EPSG:4326)
# geometry 参数指定几何数据，crs 参数指定坐标参考系统（Coordinate Reference System）
gdf = gpd.GeoDataFrame(geometry=[point, line, polygon], crs="EPSG:4326")

# 添加属性数据
# 创建一个新列 'name'，分别给点、线、多边形对象赋予名称
gdf['name'] = ['Point', 'Line', 'Polygon']

# 输出初始GeoDataFrame及其坐标系
# 使用 print 函数输出 GeoDataFrame 的内容和初始坐标系信息
print("Initial GeoDataFrame:")
print(gdf)
print("\nInitial CRS:", gdf.crs)

# 将GeoDataFrame投影到新的坐标系 (Mercator, EPSG:3857)
# to_crs 方法用于将 GeoDataFrame 的坐标系从一个 CRS 转换到另一个 CRS
# EPSG:4326 是 WGS84 地理坐标系，EPSG:3857 是 Web Mercator 投影坐标系
gdf_projected = gdf.to_crs("EPSG:3857")

# 输出投影后的GeoDataFrame及其坐标系
# 使用 print 函数输出投影后的 GeoDataFrame 的内容和新坐标系信息
print("\nProjected GeoDataFrame:")
print(gdf_projected) #输出投影后的GeoDataFrame的内容。
print("\nProjected CRS:", gdf_projected.crs) # 输出投影后的坐标参考系统 (CRS) 信息

四、基本用法

下面代码的主要目的是创建一个包含点、线和多边形几何对象的GeoDataFrame，并将其保存为Shapefile文件，以便进行地理空间数据的存储和后续处理。创建数据：

import geopandas as gpd  # 导入GeoPandas库，用于处理地理空间数据
from shapely.geometry import Point, LineString, Polygon  # 导入Shapely库中的几何对象类，用于创建几何对象

# 创建点几何对象，每个点由一对坐标 (x, y) 表示
points = [Point(1, 1), Point(2, 2), Point(3, 3)]

# 创建线几何对象，每条线由一组坐标点 (x, y) 组成，表示线的节点
lines = [
    LineString([(2, 0), (2, 4), (3, 4)]),  # 第一条线，有三个节点
    LineString([(0, 0), (1, 1), (1, 2)])   # 第二条线，有三个节点
]

# 创建多边形几何对象，每个多边形由一组坐标点 (x, y) 组成，表示多边形的顶点
polygons = [
    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),  # 第一个多边形，有四个顶点
    Polygon([(2, 2), (3, 2), (3, 3), (2, 3)])   # 第二个多边形，有四个顶点
]

# 创建包含点几何对象的GeoDataFrame，并指定初始坐标系 (WGS84, EPSG:4326)
# GeoDataFrame是一个特殊的DataFrame，专门用于处理地理空间数据
gdf_points = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326")
# 添加属性列 'name'，为每个点几何对象命名
gdf_points['name'] = ['Point1', 'Point2', 'Point3']
# 添加属性列 'value'，为每个点几何对象赋予一个数值
gdf_points['value'] = [10, 20, 30]

# 创建包含线几何对象的GeoDataFrame，并指定初始坐标系 (WGS84, EPSG:4326)
gdf_lines = gpd.GeoDataFrame(geometry=lines, crs="EPSG:4326")
# 添加属性列 'name'，为每条线几何对象命名
gdf_lines['name'] = ['Line1', 'Line2']
# 添加属性列 'value'，为每条线几何对象赋予一个数值
gdf_lines['value'] = [100, 200]

# 创建包含多边形几何对象的GeoDataFrame，并指定初始坐标系 (WGS84, EPSG:4326)
gdf_polygons = gpd.GeoDataFrame(geometry=polygons, crs="EPSG:4326")
# 添加属性列 'name'，为每个多边形几何对象命名
gdf_polygons['name'] = ['Polygon1', 'Polygon2']
# 添加属性列 'value'，为每个多边形几何对象赋予一个数值
gdf_polygons['value'] = [1000, 2000]

# 定义保存Shapefile文件的路径
shapefile_path_points = 'points.shp'  # 保存点几何对象的Shapefile路径
shapefile_path_lines = 'lines.shp'    # 保存线几何对象的Shapefile路径
shapefile_path_polygons = 'polygons.shp'  # 保存多边形几何对象的Shapefile路径

# 将包含点几何对象的GeoDataFrame保存为Shapefile文件
gdf_points.to_file(shapefile_path_points)
# 将包含线几何对象的GeoDataFrame保存为Shapefile文件
gdf_lines.to_file(shapefile_path_lines)
# 将包含多边形几何对象的GeoDataFrame保存为Shapefile文件
gdf_polygons.to_file(shapefile_path_polygons)

# 输出保存路径，确认文件已保存
print(f"Shapefile已保存到: {shapefile_path_points}, {shapefile_path_lines}, {shapefile_path_polygons}")

4.1.读取地理空间数据

GeoPandas可以轻松读取多种地理空间数据格式。要读取 polygons.shp 文件的内容，可以使用GeoPandas的 read_file 函数。以下是逐行添加注释的代码示例，展示了如何读取并打印 polygons.shp 文件的内容：

import geopandas as gpd  # 导入GeoPandas库，用于处理地理空间数据

# 定义Shapefile文件的路径
shapefile_path_polygons = 'polygons.shp'  # 保存多边形几何对象的Shapefile路径

# 读取Shapefile文件，创建GeoDataFrame
gdf_polygons = gpd.read_file(shapefile_path_polygons)

# 输出GeoDataFrame的内容
print("读取的GeoDataFrame内容:")
print(gdf_polygons)

# 输出GeoDataFrame的坐标系
print("\nGeoDataFrame的坐标系:")
print(gdf_polygons.crs)

4.2.空间操作

GeoPandas提供了丰富的空间操作功能，可以执行各种地理空间分析。以下是一些示例操作：

  • 计算两个地理空间对象的交集：
intersection = gdf1.intersection(gdf2)

计算地理空间

  • 对象的缓冲区：
buffered = gdf.buffer(distance=100)
  • 计算地理空间对象之间的距离：
distance = gdf1.distance(gdf2)
  • 判断一个地理空间对象是否包含另一个对象：
contains = gdf1.contains(gdf2)

案例代码如下：

import geopandas as gpd  # 导入GeoPandas库，用于处理地理空间数据
from shapely.geometry import Point, Polygon  # 导入Shapely库中的几何对象类，用于创建几何对象

"""
创建示例几何对象
Polygon([(x1, y1), (x2, y2), ...])：创建多边形几何对象，由一组顶点坐标点组成。
Point(x, y)：创建点几何对象，由一对坐标 (x, y) 表示。
"""
polygon1 = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])  # 第一个多边形
polygon2 = Polygon([(1, 1), (3, 1), (3, 3), (1, 3)])  # 第二个多边形
point = Point(1, 1)  # 一个点

# 创建包含几何对象的GeoDataFrame
# gpd.GeoDataFrame(geometry=[几何对象], crs="EPSG:4326")：创建包含几何对象的GeoDataFrame，并指定初始坐标系为WGS84（EPSG:4326）。
gdf1 = gpd.GeoDataFrame(geometry=[polygon1], crs="EPSG:4326")
gdf2 = gpd.GeoDataFrame(geometry=[polygon2], crs="EPSG:4326")
gdf_point = gpd.GeoDataFrame(geometry=[point], crs="EPSG:4326")

# 计算两个地理空间对象的交集
# gdf1.overlay(gdf2, how='intersection')：计算两个GeoDataFrame的交集，并创建一个新的GeoDataFrame。
intersection = gdf1.overlay(gdf2, how='intersection')
print("两个多边形的交集:")
print(intersection)

# 计算地理空间对象的缓冲区
buffered = gdf1.buffer(distance=0.5)  # 计算第一个多边形的缓冲区，缓冲距离为0.5单位。
print("\n第一个多边形的缓冲区:")
print(buffered)

# 计算地理空间对象之间的距离
distance = gdf1.distance(gdf_point)  # gdf1.distance(gdf_point)：计算第一个多边形与点之间的距离。
print("\n第一个多边形与点之间的距离:")
print(distance)

# 判断一个地理空间对象是否包含另一个对象
contains = gdf1.contains(gdf_point)
print("\n第一个多边形是否包含点:")
print(contains)

4.3.数据可视化

GeoPandas具备数据可视化功能，可以绘制地理空间数据的地图。以下是一个示例，展示如何绘制GeoDataFrame的地图：

pip install matplotlib -i https://mirrors.aliyun.com/pypi/simple/

案例代码：

import geopandas as gpd  # 导入GeoPandas库，用于处理和可视化地理空间数据
import matplotlib.pyplot as plt  # 导入Matplotlib库，用于创建图形和绘图

# 创建示例几何对象
from shapely.geometry import Point, LineString, Polygon  # 导入Shapely库中的几何对象类，用于创建几何对象

# 创建几何对象
point = Point(1, 1)  # 创建一个点对象
line = LineString([(2, 0), (2, 4), (3, 4)])  # 创建一条线对象
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])  # 创建一个多边形对象

# 创建包含几何对象的GeoDataFrame，并指定初始坐标系 (WGS84, EPSG:4326)
gdf = gpd.GeoDataFrame(geometry=[point, line, polygon], crs="EPSG:4326")

# 添加属性数据
gdf['name'] = ['Point', 'Line', 'Polygon']
gdf['value'] = [10, 20, 30]

# 绘制GeoDataFrame的地图
gdf.plot()  # 使用GeoDataFrame的plot方法绘制地图

# 显示绘制的地图
plt.show()  # 使用Matplotlib的show方法显示图形

实现效果如下：

在这个示例中，使用plot方法绘制了GeoDataFrame的地图，并使用Matplotlib进行显示。

五、功能特性

5.1.空间连接

GeoPandas支持空间连接操作，可以将两个GeoDataFrame按照空间关系进行连接。这对于分析地理空间数据之间的关联非常有用。以下是一个示例，展示如何执行空间连接：

import geopandas as gpd  # 导入GeoPandas库，用于处理和可视化地理空间数据
from shapely.geometry import Point, Polygon  # 导入Shapely库中的几何对象类，用于创建几何对象

# 创建第一个GeoDataFrame，包含一些点
gdf_points = gpd.GeoDataFrame({
    'id': [1, 2, 3],  # 'id'列表示每个点的唯一标识符
    'geometry': [Point(1, 1), Point(2, 2), Point(3, 3)]  # 'geometry'列包含点的几何对象
}, crs="EPSG:4326")  # 创建一个包含点的GeoDataFrame，并指定坐标系为WGS84（EPSG:4326）

# 创建第二个GeoDataFrame，包含一些多边形
gdf_polygons = gpd.GeoDataFrame({
    'id': ['A', 'B'],  # 'id'列表示每个多边形的唯一标识符
    'geometry': [Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]), Polygon([(2, 2), (4, 2), (4, 4), (2, 4)])]  # 'geometry'列包含多边形的几何对象
}, crs="EPSG:4326")  # 创建一个包含多边形的GeoDataFrame，并指定坐标系为WGS84（EPSG:4326）

# 执行空间连接操作，将点GeoDataFrame与多边形GeoDataFrame进行空间连接
joined_gdf = gpd.sjoin(gdf_points, gdf_polygons, how="inner", predicate="within")

# 'gpd.sjoin' 函数用于进行空间连接操作
# 'gdf_points' 是第一个GeoDataFrame，包含点
# 'gdf_polygons' 是第二个GeoDataFrame，包含多边形
# 'how="inner"' 表示使用内连接，只保留在两个GeoDataFrame中都存在的几何对象
# 'predicate="within"' 表示空间关系为点在多边形内部

# 输出空间连接结果
print(joined_gdf)  # 打印连接后的GeoDataFrame

5.2.空间缓冲区分析

GeoPandas可以执行空间缓冲区分析，用于分析地理空间对象的邻近关系。可以计算一个地理空间对象的缓冲区，以确定其周围的特定距离范围内有哪些其他对象。以下是一个示例，展示如何计算地理空间对象的缓冲区：

# 计算地理空间对象的缓冲区
buffered = gdf.buffer(distance=100)

完整案例代码如下：

import geopandas as gpd  # 导入GeoPandas库，用于处理和可视化地理空间数据
from shapely.geometry import Point, Polygon  # 导入Shapely库中的几何对象类，用于创建几何对象

# 创建一个GeoDataFrame，包含一些点
gdf = gpd.GeoDataFrame({
    'id': [1, 2, 3],  # 'id'列表示每个点的唯一标识符
    'geometry': [Point(1, 1), Point(2, 2), Point(3, 3)]  # 'geometry'列包含点的几何对象
}, crs="EPSG:4326")  # 创建一个包含点的GeoDataFrame，并指定坐标系为WGS84（EPSG:4326）

# 计算每个点的缓冲区，距离为1单位
gdf['buffer'] = gdf.geometry.buffer(1)

# 'buffer' 函数用于计算几何对象的缓冲区
# '1' 表示缓冲区的距离，单位与几何对象的坐标系一致

# 输出包含缓冲区的GeoDataFrame
print(gdf)

# 绘制原始点和缓冲区
ax = gdf.plot(color='blue', edgecolor='black')  # 绘制点
gdf.set_geometry('buffer').plot(ax=ax, color='red', alpha=0.5, edgecolor='black')  # 绘制缓冲区
plt.show()  # 显示绘制的地图

效果如下：

但是会报错：

这个警告是因为在地理坐标系（经纬度）中执行缓冲区计算会导致不准确的结果。为了解决这个问题，我们需要将数据投影到平面坐标系，然后进行缓冲区计算。下面是修改后的代码，它将GeoDataFrame投影到Mercator坐标系（EPSG:3857）进行缓冲区计算，然后再将其投影回原来的地理坐标系：

import geopandas as gpd  # 导入GeoPandas库，用于处理和可视化地理空间数据
from shapely.geometry import Point  # 导入Shapely库中的几何对象类，用于创建几何对象
import matplotlib.pyplot as plt  # 导入Matplotlib库，用于绘制图形

# 创建一个GeoDataFrame，包含一些点
gdf = gpd.GeoDataFrame({
    'id': [1, 2, 3],  # 'id'列表示每个点的唯一标识符
    'geometry': [Point(1, 1), Point(2, 2), Point(3, 3)]  # 'geometry'列包含点的几何对象
}, crs="EPSG:4326")  # 创建一个包含点的GeoDataFrame，并指定坐标系为WGS84（EPSG:4326）

# 将GeoDataFrame从地理坐标系（WGS84）投影到一个平面坐标系（Mercator, EPSG:3857）
gdf_projected = gdf.to_crs("EPSG:3857")

# 计算每个点的缓冲区，距离为1000米（因为Mercator坐标系的单位是米）
gdf_projected['buffer'] = gdf_projected.geometry.buffer(1000)

# 将GeoDataFrame投影回地理坐标系（WGS84）
gdf_buffered = gdf_projected.to_crs("EPSG:4326")

# 输出包含缓冲区的GeoDataFrame
print(gdf_buffered)

# 绘制原始点和缓冲区
ax = gdf.plot(color='blue', edgecolor='black')  # 绘制原始点
gdf_buffered.set_geometry('buffer').plot(ax=ax, color='red', alpha=0.5, edgecolor='black')  # 绘制缓冲区
plt.show()  # 显示绘制的地图

如下：

5.3.数据可视化

GeoPandas内置了数据可视化功能，可以绘制各种地理空间数据的地图。可以自定义地图样式、添加标签和图例等。

示例：绘制地理空间数据的地图，并自定义样式、添加标签和图例

import geopandas as gpd  # 导入GeoPandas库，用于处理和可视化地理空间数据
from shapely.geometry import Point  # 导入Shapely库中的几何对象类，用于创建几何对象
import matplotlib.pyplot as plt  # 导入Matplotlib库，用于绘制图形
from matplotlib.patches import Patch  # 导入Patch类，用于自定义图例

# 创建一个GeoDataFrame，包含一些点
gdf = gpd.GeoDataFrame({
    'id': [1, 2, 3],  # 'id'列表示每个点的唯一标识符
    'name': ['Point 1', 'Point 2', 'Point 3'],  # 'name'列表示点的名称
    'geometry': [Point(1, 1), Point(2, 2), Point(3, 3)]  # 'geometry'列包含点的几何对象
}, crs="EPSG:4326")  # 创建一个包含点的GeoDataFrame，并指定坐标系为WGS84（EPSG:4326）

# 将GeoDataFrame从地理坐标系（WGS84）投影到一个平面坐标系（Mercator, EPSG:3857）
gdf_projected = gdf.to_crs("EPSG:3857")

# 计算每个点的缓冲区，距离为1000米（因为Mercator坐标系的单位是米）
gdf_projected['buffer'] = gdf_projected.geometry.buffer(1000)

# 将GeoDataFrame投影回地理坐标系（WGS84）
gdf_buffered = gdf_projected.to_crs("EPSG:4326")

# 创建一个Matplotlib的子图
fig, ax = plt.subplots()

# 自定义地图样式：设置点的颜色（青色）、边界线颜色（黑色）和线宽
gdf.plot(ax=ax, color='cyan', edgecolor='black', markersize=100, linewidth=1, label='Original Points')

# 绘制缓冲区，设置缓冲区的颜色（黄色）、边界线颜色（红色）、透明度和线宽
gdf_buffered.set_geometry('buffer').plot(ax=ax, color='yellow', alpha=0.5, edgecolor='red', linewidth=2, label='Buffer Zones')

# 添加标签
for x, y, label in zip(gdf.geometry.x, gdf.geometry.y, gdf['name']):
    ax.text(x, y, label, fontsize=12, ha='right')  # 在每个点旁边添加标签，字体大小为12，水平对齐方式为右对齐

# 自定义图例处理器，设置图例项的颜色和标签
handles = [Patch(color='cyan', label='Original Points'), Patch(color='yellow', alpha=0.5, edgecolor='red', label='Buffer Zones')]

# 添加图例
plt.legend(handles=handles)

# 设置图标题
plt.title('Custom Styled Geospatial Data Visualization with Buffers')

# 显示绘制的地图
plt.show()

该代码将生成一个包含点和缓冲区的地图。点以青色显示，缓冲区以黄色半透明显示，并有红色边界线。地图中添加了每个点的标签和自定义的图例。

 

posted @ 2024-07-05 22:06  酒剑仙*  阅读(579)  评论(0)    收藏  举报