（数据科学学习手札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中的数据结构。
　　geopandas的安装和使用须要若干依赖包，若是不事先妥善安装好这些依赖包而直接使用pip install geopandas或conda install geopandas可能会引起依赖包相关错误致使安装失败，官方文档中的推荐安装方式为：

conda install --channel conda-forge 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中的投影坐标系管理，敬请期待。