1 geopandas库简介与安装
geopandas是建立在GEOS、GDAL、PROJ等开源地理空间计算相关框架之上的,类似pandas语法风格的空间数据分析Python库,其目标是尽可能地简化Python中的地理空间数据处理,减少对Arcgis、PostGIS等工具的依赖,使得处理地理空间数据变得更加高效简洁,打造纯Python式的空间数据处理工作流。
安装环境推荐是python3.8,可以使用pip方式安装:
pip install shapely==2.0.0
pip install geopandas==0.13.2
亦可以使用conda方式安装:
conda install shapely=2.0.0
conda install geopandas=0.13.2
运行完成之后,在jupyter lab中来查看geopandas的版本信息:
import geopandas as gpd
import shapely
from shapely import geometry
print("shapely版本:", shapely.__version__)
print("geopandas版本:", gpd.__version__)
2 geopandas数据结构
geopandas作为pandas向地理分析计算方面的延拓,基础的数据结构延续了Series和DataFrame的特点,创造出GeoSeries与GeoDataFrame两种基础数据结构:
2.1 geopandas.GeoSeries
2.1.1 GeoSeries中的基础几何对象
与Series相似,GeoSeries用来表示一维向量,只不过这里的向量每个位置上的元素都表示着一个shapely中的几何对象,有如下几种类型:
- Point
对应shapely.geometry中的Point,用于表示单个点,下面创建一个由若干Point对象组成的GeoSeries并像Series一样定义索引:
# 创建存放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'])
可以看到创建出的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'])
在jupyter lab中可以直接显示GeoSeries中的单个元素的矢量图:
# 直接显示GeoSeries单个元素的矢量图
gpd.GeoSeries([geometry.MultiPoint([(0, 1), (1, 0)]),
geometry.MultiPoint([(0, 0), (1, 1)])],
index=['a', 'b'])[0]
- 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'])
- 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'])
- 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'])
- 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))])])
- 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'])
- 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'])
在同一个GeoSeries可以混合上述类型中的多种几何对象,这意味着点线面在概念上相异的几何对象可以共存于同一份数据中。
2.1.2 GeoSeries常用属性
类似pandas中的Series,GeoSeries在被创建完成之后也拥有很多实用的地理属性,下面对其中较为常用的进行列举:
- is_valid
在shapely中涉及到很多拓扑计算操作时,对几何对象的合法性有要求,譬如定义多边形时坐标按顺序连线时穿过了之前定义的边就属于非法,因为geopandas对矢量对象的计算依赖于shapely,于是引进了属性用于判断每个几何对象是否合法。
下面创建两个形状相同的多边形,其中一个满足上述所说的非法情况,另一个由两个多边形拼接而成:
# 创建两个形状相同的多边形,其中一个满足上述所说的非法情况,另一个由两个多边形拼接而成
gs = 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)])])])
从形状上看两者相同,但是矢量图却显示不同的颜色:
尝试用shapely中的intersection方法来取得这两个几何对象的相交部分,出现了拓扑逻辑错误:
gs[0].intersection(gs[1])
查看gs.is_valid,可以看出第一个自相交的多边形非法:
- geom_type
# 创建混合点线面的GeoSeries,这里第5个有孔多边形内部空洞创建时使用[::-1]颠倒顺序
# 是因为GeoSeries.plot()方法绘制有孔多边形的一个bug,即外部边框与内部孔洞创建时,
# 坐标方向同为顺时针或顺时针时内部孔洞会自动被填充。
# 如果你对这个bug感兴趣,可以前往https://github.com/geopandas/geopandas/issues/951查看细节
gs = 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
gs.plot() # 对gs进行简单的可视化
geom_type返回每个几何对象类型:
-
area和length
area和length属性分别返回与GeoSeries中每个元素对应的面积值和边长值(这里的面积单位和下长度单位取决于投影坐标系,实际上,要获取正确的面积和边长,需要设置合适的投影坐标系,之后关于geopandas投影坐标系管理的文章将会详细介绍。这里仅做演示)
-
centroid和bounds
centroid属性返回每个几何对象的重心(几何中心),bounds属性返回每个几何对象所在box左下角、右上角的坐标:
-
boundary
boundary属性返回每个几何对象的低维简化表示(点对象无具体的更低维简化,故无返回值):
-
exterior与interiors
对于多边形对象,exterior返回LinearRing格式的外边框线;对于有孔多边形,interiors返回所有内部孔洞LinearRing格式边框线集合:
-
convex_hull
convex_hull属性返回每个几何对象的凸包,Polygon格式,即恰巧包含对应几何对象的凸多边形:
import numpy as np
# 利用独立的正态分布随机数创建两个MultiPoint集合
gs = 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())])
gs
ax = gs.plot(color='red') # 绘制gs
gs.convex_hull.plot(ax=ax, alpha=0.4) # 叠加绘制各自对应凸包,调低填充透明度以显示更明显
- envelope
envelope属性返回对应几何对象的box范围,Polygon格式,即包含对应元素中所有点的最小矩形:
ax = gs.plot(color='red') # 绘制gs
gs.envelope.plot(ax=ax, alpha=0.4) # 叠加绘制各自对应envelope,调低填充透明度以显示更明显
2.2 geopandas.GeoDataFrame
2.2.1 GeoDataFrame基础
顾名思义,geopandas中的GeoDataFrame是在pandas.DataFrame的基础上,加入空间分析相关内容进行改造而成,其最大特点在于其在原有数据表格基础上增加了一列GeoSeries使得其具有矢量性,所有对于GeoDataFrame施加的空间几何操作也都作用在这列指定的几何对象之上。
下面举个简单的例子,基于不同均值和标准差的正态分布随机数,创建GeoDataFrame来记录这些信息:
# geopandas.GeoDataFrame
contents = [(loc, 0.5) for loc in range(0, 10, 2)]
gdf = gpd.GeoDataFrame(data=contents,
geometry=[geometry.MultiPoint(np.random.normal(loc=loc, scale=scale, size=[10, 2]).tolist())
for loc, scale in contents],
columns=['均值', '标准差'])
gdf
其中定义GeoDataFrame时作为每行所关联几何对象的GeoSeries需要通过geometry参数指定,而除了用上述的方式创建GeoDataFrame,先创建数据表,再添加矢量信息列亦可,这时几何对象列的名称可以自由设置,但一定要利用GeoDataFrame.set_geometry() 方法将后添加的矢量列指定为矢量主列,因为每个GeoDataFrame若在定义之处没有指定矢量列,后将无法进行与适量信息挂钩的所有操作(GeoSeries所有属性都可同样作用于GeoDataFrame,因为所有空间操作实际上都直接作用于其矢量主列):
- 添加矢量列但未定义
contents = [(loc, 0.5) for loc in range(0, 10, 2)]
gdf = gpd.GeoDataFrame(data=contents, columns=['均值', '标准差'])
gdf['raw_points'] = [geometry.MultiPoint(np.random.normal(loc=loc, scale=scale, size=[10, 2]).tolist())
for loc, scale in contents]
# 尝试查看矢量类型
gdf.geom_type
这时所有直接针对GeoDataFrame的矢量相关操作都无法使用。
- 重新为GeoDataFrame指定矢量列
gdf.set_geometry('raw_points', inplace=True)
# 尝试查看矢量类型
gdf.geom_type
这时相关操作可正常使用。
- 多个矢量列切换
# 绘制第一图层
ax = gdf.plot(color='red')
gdf['convex_hull'] = gdf.convex_hull
gdf.set_geometry('convex_hull', inplace=True)
# 绘制第二图层
gdf.plot(ax=ax, color='blue', alpha=0.4)
2.2.2 GeoDataFrame数据索引
作为pandas.DataFrame的延伸,GeoDataFrame同样支持pandas.DataFrame中的.loc以及.iloc对数据在行、列尺度上进行索引和筛选,这里我们以geopandas自带的世界地图数据为例:
# 读取世界地图
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.plot()
-
查看其表格内容:
-
使用.loc和条件筛选数据:
-
使用.iloc筛选数据:
-
而除了这些常规的数据索引方式之外,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')
参考文献: