（数据科学学习手札74）基于geopandas的空间数据分析——数据结构篇

本文对应代码已上传至我的Github仓库https://github.com/CNFeffery/DataScienceStudyNotes

1 简介

　　geopandas是建立在GEOS、GDAL、PROJ等开源地理空间计算相关框架之上的，类似pandas语法风格的空间数据分析Python库，其目标是尽可能地简化Python中的地理空间数据处理，减少对Arcgis、PostGIS等工具的依赖，使得处理地理空间数据变得更加高效简洁，打造纯Python式的空间数据处理工作流。本系列文章就将围绕geopandas及其使用过程中涉及到的其他包进行系统性的介绍说明，每一篇将尽可能全面具体地介绍geopandas对应方面的知识，计划涵盖geopandas的数据结构、投影坐标系管理、文件IO、基础地图制作、集合操作、空间连接与聚合。

　　作为基于geopandas的空间数据分析系列文章的第一篇，通过本文你将会学习到geopandas中的数据结构。

　　推荐在3.9版本的Python环境中使用geopandas，注意，geopandas的安装和使用需要若干依赖包，如果不事先妥善安装好这些依赖包而直接使用pip install geopandas或conda install geopandas可能会引发依赖包相关错误导致安装失败，官方文档中的推荐安装方式为：

conda install --channel conda-forge geopandas

　　而鉴于广大缺乏“特殊上网技巧”的读者朋友们在安装geopandas时经常遇到的重重困难，我特在下方不定期更新最新稳定的geopandas安装方法（适用于WIN10及以上版本，因为我身边没有WIN7系统电脑，且Linux与Mac上的安装直接pip就可以无需专门撰写教程），请你在网络顺畅的前提下，直接在控制台中粘贴以下所有代码一口气运行即可（如无需开辟新的虚拟环境，则自行跳过前两行代码）：

• 最近更新： 2024-04-17

conda create -n geopandas-env python=3.9 -c https://mirrors.sustech.edu.cn/anaconda/pkgs/main/ -y
conda activate geopandas-env
conda install geopandas pyogrio -c https://mirrors.sustech.edu.cn/anaconda/cloud/conda-forge -y

　　如果你在使用conda执行上述最后一行命令时一直卡壳，那么可以换成使用mamba来执行安装，具体命令如下：

conda install mamba -c conda-forge -y

mamba install geopandas pyogrio -c https://mirrors.sustech.edu.cn/anaconda/cloud/conda-forge -y

　　这一切会运行的很快，运行完成之后，建议在当前虚拟环境下安装jupyter notebook或jupyter lab来查看geopandas的版本信息，像我这样就恭喜你配置完成啦😇：

　　conda-forge是一个社区项目，在conda的基础上提供了更广泛更丰富的软件资源包，通过它我们可以自动下载安装好所有geopandas的必要依赖包而无需手动繁琐地去安装它们。在完成安装后，下面我们开始对geopandas的系统性学习之旅。

2 数据结构

　　geopandas作为pandas向地理分析计算方面的延拓，基础的数据结构延续了Series和DataFrame的特点，创造出GeoSeries与GeoDataFrame两种基础数据结构：

2.1 GeoSeries

2.1.1 GeoSeries中的基础几何对象

　　与Series相似，GeoSeries用来表示一维向量，只不过这里的向量每个位置上的元素都表示着一个shapely中的几何对象，有如下几种类型：

• Points

　　对应shapely.geometry中的Point，用于表示单个点，下面我们创建一个由若干Point对象组成的GeoSeries并像Series一样定义索引：

from shapely import geometry
import geopandas as gpd

# 创建存放Point对象的GeoSeries
# 这里shapely.geometry.Point(x, y)用于创建单个点对象
gpd.GeoSeries([geometry.Point(0, 0),
               geometry.Point(0, 1),
               geometry.Point(1, 1),
               geometry.Point(1, 0)],
              index=['a', 'b', 'c', 'd'])

图1

　　可以看到创建出的GeoSeries数据类型为geometry，即几何对象。

• MultiPoint

　　对应shapely中的MultiPoint，用于表示多个点的集合，下面我们创建一个由若干MultiPoint对象组成的GeoSeries：

# 创建存放MultiPoint对象的GeoSeries
# 这里shapely.geometry.MultiPoint([(x1, y1), (x2, y2), ...])用于创建多点集合
gpd.GeoSeries([geometry.MultiPoint([(0, 1), (1, 0)]),
               geometry.MultiPoint([(0, 0), (1, 1)])],
              index=['a', 'b'])

图2

　　在jupyter notebook或jupyter lab中可以图像的形式直接显示GeoSeries中的单个元素：

图3

• LineString

　　对应shapely中的LineString，用于表示由多个点按顺序连接而成的线，下面我们创建一个由若干LineString对象组成的GeoSeries：

# 创建存放LineString对象的GeoSeries
# 这里shapely.geometry.LineString([(x1, y1), (x2, y2), ...])用于创建多点按顺序连接而成的线段
gpd.GeoSeries([geometry.LineString([(0, 0), (1, 1), (1, 0)]),
               geometry.LineString([(0, 0), (0, 1), (-1, 0)])],
              index=['a', 'b'])

图4

　　同样地，直接显示第一个元素：

图5

• MultiLineString

　　对应shapely中的MultiLineString，用于表示多条线段的集合，下面我们创建一个由若干MultiLineString对象组成的GeoSeries：

# 创建存放MultiLineString对象的GeoSeries
# 这里shapely.geometry.MultiLineString([LineString1, LineString2])用于创建多条线段的集合
gpd.GeoSeries([geometry.MultiLineString([[(0, 0), (1, 1), (1, 0)],
                                        [(-0.5, 0), (0, 1), (-1, 0)]])],
              index=['a'])

图6

　　同样地，直接显示第一个元素：

图7

• Polygon(无孔)

　　geopandas中的Polygon对应shapely中的Polygon，用于表示面，根据内部有无孔洞可继续细分。下面我们创建一个由无孔Polygon对象组成的GeoSeries：

# 创建存放无孔Polygon对象的GeoSeries
# 这里shapely.geometry.Polygon([(x1, y1), (x2, y2),...])用于创建无孔面
gpd.GeoSeries([geometry.Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])],
              index=['a'])

图8

　　同样地，直接显示第一个元素：

图9

• Polygon(有孔)

　　区分于上文中的无孔Polygon，下面我们创建一个由有孔Polygon对象组成的GeoSeries：

# 创建存放有孔Polygon对象的GeoSeries
# 这里shapely.geometry.Polygon(polygonExteriors, interiorCoords)用于创建有孔面
# 其中polygonExteriors用于定义整个有孔Polygon的外围，是一个无孔的多边形
# interiorCoords是用于定义内部每个孔洞（本质上是独立的多边形）的序列
gpd.GeoSeries([geometry.Polygon([(0,0),(10,0),(10,10),(0,10)], 
                                [((1,3),(5,3),(5,1),(1,1)), 
                                 ((9,9),(9,8),(8,8),(8,9))])])

　　同样地，直接显示第一个元素：

图10

• MultiPolygon

　　对应shapely中的MultiPolygon，用于表示多个面的集合，下面我们创建一个由MultiPolygon对象组成的GeoSeries：

# 创建存放MultiPolygon对象的GeoSeries
# 这里shapely.geometry.MultiPolygon([Polygon1, Polygon2])用于创建多个面的集合
gpd.GeoSeries([geometry.MultiPolygon([geometry.Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]),
                                      geometry.Polygon([(2, 2), (2, 3), (3, 3), (3, 2), (2, 2)])])],
              index=['a'])

图11

　　显示第一个元素：

图12

• LinearRing

　　LinearRing对应shapely.geometry中的LinearRing，是一种特殊的几何对象，可以理解为闭合的线或无孔多边形的边框，创建时传入数据的格式与Polygon相同，下面我们创建一个由LinearRing对象组成的GeoSeries：

# 创建存放LinearRing对象的GeoSeries
# 这里shapely.geometry.LinearRing([(x1, y1), (x2, y2),...])用于创建LinearRing
gpd.GeoSeries([geometry.LinearRing([(0, 0), (0, 1), (1, 1), (1, 0)])],
              index=['a'])

图13

　　显示第一个元素，可以看出LinearRing就是无孔多边形的边框线：

图14

　　在同一个GeoSeries可以混合上述类型中的多种几何对象，这意味着点线面在概念上相异的几何对象可以共存于同一份数据中

2.1.2 GeoSeries常用属性

　　类似pandas中的Series，GeoSeries在被创建完成之后也拥有很多实用的地理属性，下面对其中较为常用的进行列举：

• area

　　area属性返回与GeoSeries中每个元素一一对应的面积值（这里的面积单位和下文涉及的长度单位取决于投影坐标系，之后关于geopandas投影坐标系管理的文章将会详细介绍，这里仅做演示）：

# 创建混合点线面的GeoSeries，这里第5个有孔多边形内部空洞创建时使用[::-1]颠倒顺序
# 是因为GeoSeries.plot()方法绘制有孔多边形的一个bug，即外部边框与内部孔洞创建时坐标
# 方向同为顺时针或顺时针时内部孔洞会自动被填充，如果你对这个bug感兴趣，可以前往
# https://github.com/geopandas/geopandas/issues/951查看细节
s = gpd.GeoSeries([geometry.Polygon([(0, 0), (0.5, 0.5), (1, 0), (0.5, -0.5)]),
                   geometry.Polygon([(1, 1), (1.5, 1.5), (2, 1), (1.5, -1.5)]),
                   geometry.Point(3, 3),
                   geometry.LineString([(2, 2), (0, 3)]),
                   geometry.Polygon([(4, 4), (8, 4), (8, 8), (4, 8)], 
                                [[(5, 5), (7, 5), (7, 7), (5, 7)][::-1]])])

# 在jupyter中开启matplotlib交互式绘图模式
%matplotlib widget
s.plot() # 对s进行简单的可视化

图15

　　可以看到，s中包含了多种几何对象，下面直接得到s的面积：

图16 计算GeoSeries面积

• bounds

　　bounds属性返回每个几何对象所在box左下角、右上角的坐标信息：

图17

• length

　　length属性返回每个几何对象边长：

图18

• geom_type

　　geom_type返回每个几何对象类型：

图19

• exterior与interiors

　　对于多边形对象，exterior返回LinearRing格式的外边框线，对于有孔多边形，interiors返回所有内部孔洞LinearRing格式边框线集合：

图20

• is_valid

　　在shapely中涉及到很多拓扑计算操作时，对几何对象的合法性有要求，譬如定义多边形时坐标按顺序连线时穿过了之前定义的边就属于非法，因为geopandas对矢量对象的计算依赖于shapely，于是引进了属性用于判断每个几何对象是否合法，下面我们创建两个形状相同的多边形，其中一个满足上述所说的非法情况，另一个由两个多边形拼接而成：

s_ = gpd.GeoSeries([geometry.Polygon([(4, 0), (6, 1), (4, 1), (6, 0)]),
               geometry.MultiPolygon([geometry.Polygon([(4, 0), (5, 0.5), (6, 0)]), 
                                      geometry.Polygon([(5, 0.5), (6, 1), (4, 1)])])])

　　从形状上看两者相同：

图21

　　下面我们尝试用shapely中的intersection方法来取得这两个几何对象的相交部分，出现了拓扑逻辑错误：

图22

　　查看s_.is_valid，可以看出第一个自相交的多边形非法：

图23

• boundary

　　boundary返回每个几何对象的低维简化表示（点对象无具体的更低维简化，故无返回值）：

图24

• centroid

　　centroid返回每个几何对象的重心（几何中心）：

图25

• convex_hull

　　convex_hull返回每个几何对象的凸包，Polygon格式，即恰巧包含对应几何对象的凸多边形：

import numpy as np

# 利用独立的正态分布随机数创建两个MultiPoint集合
s__ = gpd.GeoSeries([geometry.MultiPoint(np.random.normal(loc=0, scale=2, size=[10, 2]).tolist()),
               geometry.MultiPoint(np.random.normal(loc=5, scale=2, size=[10, 2]).tolist())])

ax = s__.plot(color='red') # 绘制s__
s__.convex_hull.plot(ax=ax, alpha=0.4) # 叠加绘制各自对应凸包，调低填充透明度以显示更明显

图26

• envelope

　　envelope属性返回对应几何对象的box范围，Polygon格式，即包含对应元素中所有点的最小矩形：

import numpy as np

# 创建两团独立的MultiPoint
s__ = gpd.GeoSeries([geometry.MultiPoint(np.random.normal(loc=0, scale=2, size=[10, 2]).tolist()),
               geometry.MultiPoint(np.random.normal(loc=5, scale=2, size=[10, 2]).tolist())])

ax = s__.plot(color='red') # 绘制s__
s__.envelope.plot(ax=ax, alpha=0.4) # 叠加绘制各自对应envelope，调低填充透明度以显示更明显

图27

2.2 GeoDataFrame

2.2.1 GeoDataFrame基础

　　顾名思义，geopandas中的GeoDataFrame是在pandas.DataFrame的基础上，加入空间分析相关内容进行改造而成，其最大特点在于其在原有数据表格基础上增加了一列GeoSeries使得其具有矢量性，所有对于GeoDataFrame施加的空间几何操作也都作用在这列指定的几何对象之上。下面我们举个简单的例子，基于不同均值和标准差的正态分布随机数，创建GeoDataFrame来记录这些信息：

contents = [(loc, 0.5) for loc in range(0, 10, 2)]
geo_df = gpd.GeoDataFrame(data=contents, 
                          geometry=[geometry.MultiPoint(np.random.normal(loc=loc, scale=scale, size=[10, 2]).tolist()) 
                                    for loc, scale in contents],
                          columns=['均值', '标准差'])
geo_df

图28

　　其中定义GeoDataFrame时作为每行所关联几何对象的GeoSeries需要通过geometry参数指定，而除了用上述的方式创建GeoDataFrame，先创建数据表，再添加矢量信息列亦可，这时几何对象列的名称可以自由设置，但一定要利用GeoDataFrame.set_geometry()方法将后添加的矢量列指定为矢量主列，因为每个GeoDataFrame若在定义之处没有指定矢量列，后将无法进行与适量信息挂钩的所有操作（GeoSeries所有属性都可同样作用于GeoDataFrame，因为所有空间操作实际上都直接作用于其矢量主列）：

• 添加矢量列但未定义

geo_df = gpd.GeoDataFrame(contents, columns=['均值', '标准差'])
geo_df['raw_points'] = [geometry.MultiPoint(np.random.normal(loc=loc, scale=scale, size=[10, 2]).tolist()) 
                                    for loc, scale in contents]
# 尝试查看矢量类型
geo_df.geom_type

图29

　　这时所有直接针对GeoDataFrame的矢量相关操作都无法使用。

• 重新为GeoDataFrame指定矢量列

geo_df.set_geometry('raw_points').geom_type

　　这时相关操作可正常使用：

图30

• 多个矢量列切换

　　通过前面的内容，我们知道了每个GeoDataFrame都有一个矢量主列，相关操作例如绘图都基于此列，实际上GeoDataFrame允许表中存在多个矢量列，只要求任意时刻有且仅有1列为矢量主列即可，因此我们可以在一个GeoDataFrame中保存多列矢量，需要用到哪列时再进行切换即可，如下面的例子：

geo_df = gpd.GeoDataFrame(contents, columns=['均值', '标准差'])
geo_df['raw_points'] = [geometry.MultiPoint(np.random.normal(loc=loc, scale=scale, size=[10, 2]).tolist()) 
                                    for loc, scale in contents]
geo_df.set_geometry('raw_points', inplace=True) # inplace=True表示对原数据进行更新

# 绘制第一图层
ax = geo_df.plot(color='red')
geo_df['convex_hull'] = geo_df.convex_hull
# 切换矢量主列
geo_df.set_geometry('convex_hull', inplace=True)

# 绘制第二图层
geo_df.plot(ax=ax, color='blue', alpha=0.4)

图31

2.2.2 GeoDataFrame数据索引

　　作为pandas.DataFrame的延伸，GeoDataFrame同样支持pandas.DataFrame中的.loc以及.iloc对数据在行、列尺度上进行索引和筛选，这里我们以geopandas自带的世界地图数据为例：

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.plot()

图32 geopandas自带世界地图

　　查看其表格内容：

图33

　　使用.loc+条件筛选选择数据：

图34

　　使用.iloc选择数据：

图35

　　而除了这些常规的数据索引方式之外，geopandas为GeoDataFrame添加了.cx索引方式，可以传入所需的空间范围，用于索引与传入范围相交的对应数据：

# 选择与东经80度-110度，北纬0度-30度范围相交的几何对象
part_world = world.cx[80:110, 0:30]

# 绘制第一图层：世界地图
ax = world.plot(alpha=0.05)
# 绘制第二图层：.cx所选择的地区
ax = part_world.plot(ax=ax, alpha=0.6)
# 绘制第三图层：.cx条件示意图
ax = gpd.GeoSeries([geometry.box(minx=80, miny=0, maxx=110, maxy=30).boundary])\
    .plot(ax=ax, color='red')

　　示意图如下：

图36

　　放大到所选区域，可以看出正如前面所说，通过.cx，所有与指定空间范围有重叠的对象都被选择：

图37

　　以上就是本文的全部内容，如有笔误望指出，系列文章下一篇将详细介绍geopandas中的投影坐标系管理，敬请期待。