文章首页 Web开发 Windows开发 编程语言 数据库技术 移动平台 系统服务 微信 设计 布布扣 其他

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

时间：2020-02-08 21:40:13   收藏：0   阅读：149

本文对应代码已上传至我的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中的几何对象，有如下几种类型：

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对象的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对象的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对象的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对象的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对象的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对象的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对象的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在被创建完成之后也拥有很多实用的地理属性，下面对其中较为常用的进行列举：

# 创建混合点线面的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面积

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

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
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的矢量相关操作都无法使用。

geo_df.set_geometry('raw_points').geom_type

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

技术分享图片
图30

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

原文：https://www.cnblogs.com/feffery/p/11898190.html

评论（0）
© 2014 bubuko.com 版权所有 - 联系我们:wmxa8@hotmail.com
打开技术之扣，分享程序人生！