本文主要介绍使用GeoPandas的要点。 GeoPandas 是一个 Python 开源项目，旨在为地理空间数据处理提供丰富而简单的接口。 GeoPandas 扩展了 Pandas 的数据类型并使用 matplotlib 进行绘图。 GeoPandas的官方仓库地址是：GeoPandas。 GeoPandas的官方文档地址是：GeoPandas-doc。本文主要参考GeoPandas Examples Gallery。 GeoPandas的基本使用请参见Python绘制数据地图1-GeoPandas入门指南。关于GeoPandas可视化的介绍，请参见Python绘制数据地图2-GeoPandas地图可视化。

本文所有代码请参见：Python学习笔记

GeoPandas推荐使用Python 3.7及以上版本，最佳运行环境是Linux系统。 GeoPandas安装命令如下：

pip install geopandas

如果以上命令安装出现问题，建议使用conda安装GeoPandas。命令如下：

conda 安装geopandas

或：

conda install --channel conda-forge geopandas

除了GeoPandas之外，还需要安装以下第三方库：

pip install mapclassify
pip 安装 matplotlib_scalebar
pip 安装 rtree
pip 上下文安装

# 移除jupyter笔记本环境的警告
进口警告
warnings.filterwarnings("忽略")
# 检查geopandas版本
将 geopandas 导入为 gpd
gpd.__版本__

'0.13.2'

1 分层图表 Choropleth

Choropleth 是表示地理区域内数据分布的可视化图表。它将地图划分为不同的区域，并使用不同程度的颜色或阴影来显示该区域的数据值。通常，分层图表用于显示人口统计和自然资源分布等数据。分层统计图表可以帮助观察者更容易地了解数据在地理空间中的分布和变化趋势，有助于决策和规划相关工作。

导入 geopandas 作为 gpd
从 geopandas 导入 read_file

# pip install mapclassify
导入地图分类
映射分类.__版本__

'2.5.0'

#读取四川地图数据，数据来自DataV.GeoAtlas，并投影到EPSG:4573
data = www.gyballet.com_file('https://www.gyballet.com/areas_v3/bound/510000_full.json').to_crs('EPSG:4573')
数据.head()

简单的分层统计

以下代码通过方案分类统计了四川省地级市所包含的区县数量。

ax = 数据.plot(
列=“childrenNum”，
schema="QUANTILES", # 设置分层着色标准
边缘颜色='浅灰色',
k=7, # 等级数
cmap=“布鲁斯”，
图例=真实，
# 通过fmt设置位数
legend_kwds={"loc": "中左", "bbox_to_anchor": (1, 0.5),"fmt": "{:.2f}"}
）
# 显示地级市包含的区县数量
对于 data.index 中的索引：
x = data.iloc[索引].geometry.centroid.x
y = data.iloc[索引].geometry.centroid.y
名称 = data.iloc[索引]["childrenNum"]
ax.text(x, y, 名称, ha=“中心”, va=“中心”,color='红色')

#查看标签，这是分级区间
labels = [t.get_text() for t in ax.get_legend().get_texts()]
标签

[' 3.00, 5.00',
'5.00，6.00'，
'6.00，6.57'，
'6.57，7.43'，
'7.43，9.29'，
'9.29, 13.57',
'13.57, 20.00']

#检查每个分级标准和区间数，一般是左开右闭
res = mapclassify.Quantiles(data["childrenNum"], k=7)
资源

分位数
间隔计数
-----------------------
[3.00, 5.00] | 5
(5.00, 6.00] | 4
(6.00, 6.57] | 0
(6.57, 7.43] | 3(7.43, 9.29] | 3
(9.29, 13.57] | 3
(13.57, 20.00] | 3

间隔显示

ax = 数据.plot(
列=“childrenNum”，
方案=“箱线图”，
边缘颜色='k',
cmap="OrRd", # 设置分层着色标准
图例=真实，
# 使用interval设置是否显示interval间隔
legend_kwds={"loc": "中左", "bbox_to_anchor": (1, 0.5), "间隔": True}
）
# 显示地级市包含的区县数量
对于 data.index 中的索引：
x = data.iloc[索引].geometry.centroid.x
y = data.iloc[索引].geometry.centroid.y
名称 = data.iloc[索引]["childrenNum"]
ax.text(x, y, 名称, ha=“中心”, va=“中心”,color='红色')

类别展示

数据以数值分类的形式显示。下辖20个区县的地级市是成都市。

ax = 数据.plot(
列=“childrenNum”，
categorical=True, # 以数值分类的形式显示
图例=真实，
cmap =“tab20”，
# 对于分类数据，fmt设置没有用
legend_kwds={"loc": "中左", "bbox_to_anchor": (1, 0.5), "fmt": "{:.0f}"},
）
# 显示地级市包含的区县数量
对于 data.index 中的索引：
x = data.iloc[索引].geometry.centroid.x
y = data.iloc[索引].geometry.centroid.y
名称 = data.iloc[索引]["childrenNum"]
ax.text(x, y, 名称, ha=“中心”, va=“中心”,color='红色')

自定义评级

#自定义分级标准
def 自定义（值）：
# 设置ABC的三个级别
级别=无如果值 > 15：
级别 = 'A'
elif 值 > 8：
级别 = 'B'
别的：
级别 = 'C'
回报水平
#根据自定义函数映射到新列
data['level'] = data['childrenNum'].apply(custom)
数据.head()

ax = 数据.plot(
列=“级别”，
categorical=True, # 以数值分类的形式显示
图例=真实，
cmap =“冷暖”，
# 对于分类数据，fmt设置没有用
legend_kwds={"loc": "中左", "bbox_to_anchor": (1, 0.5), "fmt": "{:.0f}"},
）
# 显示地级市包含的区县数量
对于 data.index 中的索引：
x = data.iloc[索引].geometry.centroid.x
y = data.iloc[索引].geometry.centroid.y
名称 = data.iloc[索引]["childrenNum"]
ax.text(x, y, 名称, ha=“中心”, va=“中心”,color='红色')

2 从 DataFrame 创建 GeoDataFrame

基于纬度和经度数据

GeoDataFrame 有一个几何列，我们可以根据纬度和经度数据 Latitude 和 Longitude 创建它。

导入 pandas 作为 pd
# 生成有关南美城市的数据框数据
df = pd.DataFrame(
{
“城市”：[“布宜诺斯艾利斯”，“巴西利亚”，“圣地亚哥”，“波哥大”，“加拉加斯”]，
“国家”：[“阿根廷”、“巴西”、“智利”、“哥伦比亚”、“委内瑞拉”]，
“纬度”：[-34.58，-15.78，-33.45，4.60，10.48]，
“经度”：[-58.66，-47.91，-70.66，-74.08，-66.86]，
}
）
df

# 将dataframe转换为GeoDataFrame
将 geopandas 导入为 gpd
gdf = gpd.GeoDataFrame(
df，几何= gpd.points_from_xy（df.经度，df.纬度），crs =“EPSG：4326”
）
广东省

#在南美地图上展示世界 = www.gyballet.com_file(gpd.datasets.get_path("naturalearth_lowres"))
# 目标南美洲
ax = www.gyballet.com[-90:-55, -25:15].plot(color=“白色”, edgecolor=“黑色”)
# 在斧头区域绘制地图
gdf.plot(ax=ax, 颜色=“红色”)

基于WTK格式数据

WKT（众所周知的文本）是一种用于描述地理位置的数据格式。 WTK格式的数据包含点、线、多边形等地理位置信息。 WTK格式的数据可以被许多GIS软件和地理位置分析工具读取和处理。我们可以将带有 WKT 数据的 DataFrame 转换为 GeoDataframe。

df = pd.DataFrame(
{
“城市”：[“布宜诺斯艾利斯”，“巴西利亚”，“圣地亚哥”，“波哥大”，“加拉加斯”]，
“国家”：[“阿根廷”、“巴西”、“智利”、“哥伦比亚”、“委内瑞拉”]，“坐标”：[
“点（-58.66 -34.58）”，
“点（-47.91 -15.78）”，
“点（-70.66 -33.45）”，
“点（-74.08 4.60）”，
“点（-66.86 10.48）”，
],
}
）
df

#创建新列然后数据转换
df["坐标"] = gpd.GeoSeries.from_wkt(df["坐标"])gdf = gpd.GeoDataFrame(df, 几何=“坐标”, crs=“EPSG:4326”)
打印（gdf.head（））

 城市国家坐标
0 布宜诺斯艾利斯 阿根廷 点 (-58.66000 -34.58000)
1 巴西利亚 巴西点 (-47.91000 -15.78000)
2 智利圣地亚哥点 (-70.66000 -33.45000)
3 哥伦比亚波哥大 POINT (-74.08000 4.60000)
4 加拉加斯委内瑞拉点 (-66.86000 10.48000)

# 显示在南美洲地图上
世界 = www.gyballet.com_file(gpd.datasets.get_path("naturalearth_lowres"))
# 目标南美洲
ax = www.gyballet.com[-90:-55, -25:15].plot(color=“白色”, edgecolor=“黑色”)
# 在斧头区域绘制地图
gdf.plot(ax=ax, 颜色=“红色”)

3 添加比例尺

在地理空间数据分析和可视化过程中，比例尺可以帮助我们了解地图上的距离和大小关系。基于matplotlib进行可视化时，可以使用matplotlib-scalebar库添加比例尺。

简单音阶

导入 geopandas 作为 gpd
# pip install matplotlib_scalebar 安装
从 matplotlib_scalebar.scalebar 导入 ScaleBar

# nybb 是纽约排名前五的地图
nybb = www.gyballet.com_file(gpd.datasets.get_path("nybb"))
# 北美通用坐标系，坐标单位为米
nybb = www.gyballet.com_crs(32619)
nybb.head()

如下所示，创建ScaleBar对象所需的唯一参数是dx。 dx表示输入图像每个像素代表的长度，units是dx的单位。该参数的值取决于坐标参考系的单位。在之前的nybb数据集中，已经使用了epsg:32619坐标系。坐标系以米为单位。如下图，可以看到www.gyballet.com输出结果中的Axis Info项标识了参考系的单位是米。

www.gyballet.com

名称：WGS 84 / UTM 区 19N
轴信息[笛卡尔]：
- E[东]：东距（米）
- N[北]：北距（米）
使用范围：
- 名称：西经 72° 至 66°W 之间，北半球赤道至北纬 84° 之间，陆上和近海。阿鲁巴岛。巴哈马。巴西。加拿大 - 新不伦瑞克省 (NB)；拉布拉多犬;努纳武特地区；新斯科舍省 (NS)；魁北克。哥伦比亚。多明尼加共和国。格陵兰。荷属安的列斯。波多黎各。特克斯和凯科斯群岛。美国。委内瑞拉。
- 边界：(-72.0, 0.0, -66.0, 84.0)
协调操作：
- 名称：UTM区19N
- 方法：横轴墨卡托
基准面：世界大地测量系统 1984 年整体
- 椭球：WGS 84
- 本初子午线：格林威治

在下面代码中添加了比例尺和像素尺寸，该比例尺采用的是线段式表示方式，即在地图上绘制一条线段并注明该地图上该线段所代表的实际距离。

ax = nybb.plot()
# 在地图中添加比例尺和像素尺寸
scalebar =ScaleBar(dx=1,units="m")
ax.add_artist(scalebar)

确定比例尺基准长度

如下所示，以经纬度为单位的epsg:4326坐标系，其单位尺度为度（经纬度）。

# nybb为纽约五大区地图
nybb = www.gyballet.com_file(gpd.datasets.get_path("nybb"))
nybb = www.gyballet.com_crs(4326)
nybb.plot()

www.gyballet.com

Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

可以通过计算该地图中相邻两点之间的距离长度来确定比例尺基准长度，要注意的是这两点应该位于待绘制的地图中。

from shapely.geometry.point import Point
points = gpd.GeoSeries(
[Point(-73.9, 40.7), Point(-74.9, 40.7)], crs=4326
)
# 将两点转换到以米为单位的坐标系
points = www.gyballet.com_crs(32619)
# 计算点之间的距离，距离单位为坐标系的单位
distance_meters = points[0].distance(points[1])
distance_meters

84698.53985065906

nybb = www.gyballet.com_crs(4326)
ax = nybb.plot()
ax.add_artist(ScaleBar(distance_meters,"m"))

比例尺自定义

通过更改 ScaleBar 参数能够调整比例尺的显示效果，ScaleBar具体参数如下所示。这些参数的使用可以自行尝试。

scalebar = ScaleBar(
dx, # 像素和长度之间的比例尺。例如，如果一个像素代表1毫米，则dx=0.001。
units="m", # 长度单位
dimension="si-length", # 维度
label=None, # 刻度尺标签
length_fraction=None, # 刻度尺长度占比
height_fraction=None, # 刻度尺高度占比
width_fraction=None, # 刻度尺宽度占比
location=None, # 刻度尺的位置
pad=None, # 刻度尺和边框之间的间距
border_pad=None, # 刻度尺和边框之间的边距
sep=None, # 刻度尺标签和刻度之间的距离
frameon=None, # 是否显示边框
color=None, # 刻度尺和标签的颜色
box_color=None, # 边框的颜色
box_alpha=None, # 边框的透明度
scale_loc=None, # 刻度线的位置
label_loc=None, # 刻度尺标签的位置
font_properties=None, # 标签和刻度线的字体属性
label_formatter=None, # 标签的格式化函数
scale_formatter=None, # 刻度线的格式化函数
fixed_value=None, # 固定的数值
fixed_units=None, # 固定的单位
animated=False, # 是否允许动画
rotation=None, # 刻度尺的旋转角度
bbox_to_anchor=None, # bbox的锚点
bbox_transform=None, # bbox的变换
)

此外，也可以更改一像素代表的长度单位，如ScaleBar(2, dimension="si-length", units="km")表示图中1像素代表实际si-length（国际单位制）中的2km。所支持的长度单位参数如下表所示：

dimension	units
si-length	km, m, cm, um
imperial-length	in, ft, yd, mi
si-length-reciprocal	1/m, 1/cm
angle	deg

一些比例尺参数调整的示例如下

nybb = www.gyballet.com_file(gpd.datasets.get_path("nybb")).to_crs(32619)
ax = nybb.plot()
# 改变位置和方向
scale1 = ScaleBar(
dx=1,
label="Scale 1",
location="lower left",  # 位置
label_loc="left",
scale_loc="top",  # 注释文字相对于横线方向位置
)
# 改变颜色
scale2 = ScaleBar(
dx=1,
label="Scale 2",
location="center",
color="#b32400",
box_color="yellow",
box_alpha=0.8,  # 透明度
)
# 改变文字
scale3 = ScaleBar(
dx=1,
label="Scale 3",
font_properties={
"size": "large",
},
location="lower right",  # 位置
scale_formatter=lambda value, unit: f"> {value} {unit} <",
)
# 改变长度
scale4 = ScaleBar(
dx=1,
label="Scale 4",
length_fraction=0.5, # 表示刻度线占绘图区域的50%
scale_loc="top",
label_loc="left",
border_pad=1,
pad=0.25,
)
ax.add_artist(scale1)
ax.add_artist(scale2)
ax.add_artist(scale3)
ax.add_artist(scale4)

4 图层操作与几何运算

4.1 图层叠加

在geopandas中，overlay()函数是用于将两个地理图层进行叠加分析的函数。它可以进行求交集、并集、差集和对称差集等操作。overlay()函数的基本语法如下：

geopandas.overlay(layer1, layer2, how)

其中，layer1和layer2是两个geopandas地理图层对象，how是一个字符串，指定要进行的叠加操作。how参数有以下取值：

• intersection：交集
• union：并集
• difference：差集
• symmetric_difference：对称差集
• identity

在下面的示例中展示overlay函数的使用方式。

准备数据

import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon, Point

# 画一个圆
center = Point(2, 2)  # 圆心坐标
radius = 1  # 圆的半径
circle = center.buffer(radius)
gdf1 = gpd.GeoDataFrame({'geometry': circle, 'circle':[0]})
gdf1.plot()

gdf1

# 画两个正方形
square = gpd.GeoSeries([Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
Polygon([(2, 2), (4, 2), (4, 4), (2, 4)])])
gdf2 = gpd.GeoDataFrame({'geometry': square, 'square':[0,1]})
gdf2.plot()

gdf2

# 展示共同绘图结果
ax = gdf1.plot()
gdf2.plot(ax=ax)

交集intersection

# 需要pip install rtree
gdf = gpd.overlay(gdf1, gdf2, how='intersection')
gdf

gdf.plot(cmap="tab10")

并集union

# 需要pip install rtree
gdf = gpd.overlay(gdf1, gdf2, how='union')
gdf

gdf.plot(cmap="tab10")

差集difference

# 需要pip install rtree
# 提取在gdf1中，但不在gdf2中的区域
gdf = gpd.overlay(gdf1, gdf2, how='difference')
# 也可以用以下写法更加直观
# gdf = gdf1.overlay(gdf2, how='difference')
gdf

gdf.plot(cmap="tab10")

对称差集symmetric_difference

# 需要pip install rtree
# 提取不在gdf1和pdf2交集的区域
gdf = gpd.overlay(gdf1, gdf2, how='symmetric_difference')
gdf

gdf.plot(cmap="tab10")

identity

identity是ArcGIS中常用的操作。意思是将源地理图层与参考图层进行比较，以在源图层中标识与参考图层中相交的区域。使用identity的一个典型场景是当需要分析两个图层交集的时候。例如，可能有一个图层包含了所有的道路，另一个图层包含了所有的建筑。通过使用identity可以找到所有的建筑物位于哪些道路上。

# 需要pip install rtree
gdf = gpd.overlay(gdf1, gdf2, how='identity')
gdf

gdf.plot(cmap="tab10")

4.2 空间连接

空间连接允许将两个或多个空间数据集合并成一个新的数据集。例如，我们有两个数据集，一个包含所有城市的边界，另一个包含所有的人口数据。通过空间连接，我们可以将这两个数据集合并成一个新的数据集，其中每个城市都会有相应的人口数据。GeoPandas提供sjoin函数将两个GeoDataFrame数据集基于空间关系进行连接。sjoin函数常用参数如下：

sjoin(left_df, right_df, how='inner', op='intersects', lsuffix='left', rsuffix='right')

其中，参数含义如下：

• left_df：左侧的GeoDataFrame数据集。
• right_df：右侧的GeoDataFrame数据集。
• how：连接方式，可选项如下：
  • inner （默认选项）：返回两个GeoDataFrame中具有共同空间索引的几何体的交集。
  • left：返回左侧GeoDataFrame中的所有几何体，以及右侧GeoDataFrame中与之相交的几何体。如果右侧GeoDataFrame中没有与左侧相交的几何体，则右侧数据中的所有列都将为null。
  • right：与left相反，返回右侧GeoDataFrame中的所有几何体，以及左侧GeoDataFrame中与之相交的几何体。如果左侧GeoDataFrame中没有与右侧相交的几何体，则左侧数据中的所有列都将为null。
• predicate：连接的空间关系，常用选项如下：
  • intersects （默认选项）：返回两个几何体相交的所有几何体。
  • contains：返回左侧GeoDataFrame中包含于右侧GeoDataFrame中的所有几何体。
  • within：返回右侧GeoDataFrame中包含于左侧GeoDataFrame中的所有几何体。
  • touches：返回两个几何体相切的所有几何体。
  • crosses：返回两个几何体相交但不相切的所有几何体。
  • overlaps：返回两个几何体部分重叠的所有几何体。
• lsuffix：组合后左侧数据集中几何对象列的后缀，默认为left。
• rsuffix：组合后右侧数据集中几何对象列的后缀，默认为right。

以下示例展示了如何使用sjoin函数进行空间连接。

准备数据

# 创建点 GeoDataFrame
points = gpd.GeoDataFrame(
[
{'id': 'p1', 'geometry': Point(0, 0)},
{'id': 'p2', 'geometry': Point(1, 1)},
{'id': 'p3', 'geometry': Point(2, 2)},
{'id': 'p4', 'geometry': Point(3, 3)}
],
crs='EPSG:4326'
)
points

points.plot()

# 创建多边形 GeoDataFrame
polygons = gpd.GeoDataFrame(
[
{'id': 'P1', 'geometry': Polygon([(0, 0), (0, 2), (2, 2), (2, 0)])},
{'id': 'P2', 'geometry': Polygon([(1, 1), (1, 3), (3, 3), (3, 1)])}
],
crs='EPSG:4326'
)
polygons

polygons.plot()

sjoin函数

# 左连接
join_left_df = points.sjoin(polygons, how="left")
# 输出结果中的每一行都表示左侧GeoDataFrame中的一个几何对象与右侧GeoDataFrame中的一个几何对象进行了连接后得到的结果。
# index_right表示右侧GeoDataFrame中的行索引
# id_left：表示左侧GeoDataFrame中的几何对象的ID
# id_right：表示右侧GeoDataFrame中的几何对象的ID
# geometry：表示连接后的几何对象
join_left_df

# 右连接
join_right_df = points.sjoin(polygons, how="right")
join_right_df

# 设置predicate
join_right_within_df = points.sjoin(polygons, how="left", predicate="contains")
join_right_within_df

4.3 几何操作

GeoPandas提供了多种用于几何操作的函数，具体如下：

• 构造方法
  • buffer(distance, resolution=16)：返回一个GeoSeries，其中包含与每个几何对象距离在给定距离内的所有点的几何形状。
  • boundary：返回一个GeoSeries，其中包含每个几何形状的集合理论边界的低维对象。
  • centroid：返回一个GeoSeries，其中包含每个几何质心的点。
  • convex_hull：返回一个GeoSeries，其中包含表示包含每个对象中所有点的最小凸多边形的几何形状，除非对象中的点数小于三个。对于两个点，凸包会折叠成一个线串；对于一个点，凸包是一个点。
  • envelope：返回一个GeoSeries，其中包含包含每个对象的点或最小矩形多边形（其边与坐标轴平行）的几何形状。
  • simplify(tolerance, preserve_topology=True)：返回一个GeoSeries，其中包含每个对象的简化表示。在geopandas中，simplify函数可以用来简化多边形的形状，以减少地图数据的大小，同时也可以提高绘图的效率。当绘图数据特别大时，该函数很有用。tolerance：简化容差值，代表简化几何对象的形状后的最大允许误差。当 tolerance 值越小时，简化后的几何对象的形状越接近原始几何对象的形状。preserve_topology：是否保持拓扑结构，默认值为True，表示保持拓扑结构。
  • unary_union：返回一个几何形状，其中包含GeoSeries中所有几何形状的联合。
• 几何变化方法
  • affine_transform(self, matrix)：使用仿射变换矩阵来变换 GeoSeries 的几何形状。matrix 为一个包含6、12个元素的列表或元组（2d情况、3d情况）的仿射变换矩阵。关于 matrix 参数的使用需要有仿射变换的知识。
  • rotate(ngle, origin='center', use_radians=False)：旋转 GeoSeries 的坐标系。
  • scale(xfact=1.0, yfact=1.0, zfact=1.0, origin='center')：沿着(x, y, z)三个方向缩放 GeoSeries 的几何形状。
  • skew(xs=0.0, ys=0.0, origin='center', use_radians=False)：基于原点origin，沿着 x 和 y 两个方向倾斜/扭曲 GeoSeries 的几何形状。
  • translate(xoff=0.0, yoff=0.0, zoff=0.0)：平移 GeoSeries 的坐标系。

构造方法使用示例

import geopandas as gpd
# 加载数据集
world = www.gyballet.com_file(gpd.datasets.get_path('naturalearth_lowres'))
# 展示结果
world.plot()

# buffer函数
buffered = world.geometry.buffer(distance=5)
# 显示结果
buffered.plot()

# 获取几何形状边界
boundary = world.geometry.boundary
# 显示结果
boundary.plot()

# 获取几何质心
centroids = world.geometry.centroid
# 显示结果
centroids.plot(marker='*', color='green', markersize=5)

# 获取几何形状的凸包
convex_hulls = world.geometry.convex_hull
# 显示结果
convex_hulls.plot(alpha=0.5, edgecolor='k')

# 获取几何形状的外接矩形
envelopes = world.geometry.envelope
# 显示结果
envelopes.plot(alpha=0.5, edgecolor='k')

# 对几何对象进行简化处理
simplified = world.geometry.simplify(tolerance=0.1)
# 显示结果
simplified.plot(alpha=0.5, edgecolor='k')

merged = world.geometry.unary_union
# 将合并后的几何对象转换为GeoDataFrame
gdf_merged = gpd.GeoDataFrame(geometry=[merged])
# 打印后只有一行
print(gdf_merged)
gdf_merged.plot()

                                            geometry
0  MULTIPOLYGON (((-162.440 -79.281, -163.027 -78...

几何变化方法使用示例

# 读取数据集
import geopandas as gpd
nybb = www.gyballet.com_file(gpd.datasets.get_path('nybb'))
ax = nybb.plot()
# 启用科学计数法
ax.ticklabel_format(style='sci', axis='both', scilimits=(0,0))

from shapely.affinity import affine_transform
# 仿射变换
# 定义仿射变换参数
a, b, d, e, xoff, yoff = 1.5, 0.5, 0, 0.5, 1.5, 0
tmp = nybb.copy()
# 对nybb数据集中的几何对象进行仿射变换
tmp['geometry'] = tmp['geometry'].apply(lambda x: affine_transform(x, [a, b, d, e, xoff, yoff]))
# 显示变换后的nybb数据集
tmp.plot()

# 旋转
nybb_rotate = nybb.geometry.rotate(angle=45)
nybb_rotate.plot()

# 缩放
nybb_scale = nybb.geometry.scale(xfact=2, yfact=2, zfact=1)
nybb_scale.plot()

# 倾斜/扭曲
nybb_skew = nybb.geometry.skew(xs=0.1, ys=0.2, use_radians=True)
ax = nybb_skew.plot()
# 启用科学计数法
ax.ticklabel_format(style='sci', axis='both', scilimits=(0,0))

# 平移
nybb_translated = nybb.geometry.translate(xoff=100000, yoff=100000, zoff=0)
ax = nybb_translated.plot()
# 启用科学计数法
ax.ticklabel_format(style='sci', axis='both', scilimits=(0,0))

4.4 汇总

在geopandas中，dissolve函数可以对具有相同属性值的几何对象进行合并，从而生成新的几何对象。在汇总过程中，可以选择保留某些字段的信息，也可以对其他字段进行统计计算。dissolve函数如下：

geopandas.GeoDataFrame.dissolve(by=None, aggfunc='first', as_index=True, **kwargs)

函数参数介绍：

• by: 可以是一个字段名，也可以是一列字段名的列表。表示按照哪些字段进行汇总。默认为None，即将所有要素合并成一个要素。
• aggfunc: 统计函数，用于对其他字段进行计算，可以是以下函数之一：
  • 'first': 返回第一个非空值。
  • 'last': 返回最后一个非空值。
  • 'mean': 返回平均值。
  • 'sum': 返回总和。
  • 'min': 返回最小值。
  • 'max': 返回最大值。
  • 自定义函数：可以传入自定义的聚合函数。
• as_index: 是否将by参数指定的字段作为行索引，默认为True。
• *kwargs: 其他参数。

下面示例代码演示了dissolve函数的使用。

import geopandas as gpd
# 读取湖北省地图数据
data = www.gyballet.com_file('https://www.gyballet.com/areas_v3/bound/420000_full.json')
data.head()

data.plot(cmap='tab20')

# 使用dissolve函数合并几何体，根据地级市的区县数分组
dissolve_data = data.dissolve(by='childrenNum')
dissolve_data.head()

dissolve_data.plot(cmap='tab20')

# 使用dissolve函数合并几何体，根据地级市的区县数分组，其他列求均值
dissolve_data = data.dissolve(by='childrenNum', aggfunc='mean')
dissolve_data.head()

4.5 缺失值与空值处理

在使用geopandas处理地理空间数据时，经常会遇到None和Empty这两个概念。虽然它们都表示缺失值，但它们之间有着一些区别。

• None：表示属性或者列的值不存在，或者没有被填充。在geopandas中，如果一个geometry列的值为None，那意味着这个几何对象不存在。
• Empty：表示属性或者列的值存在，但是值为空。在geopandas中，如果一个geometry列的值为空，那意味着这个几何对象是存在的，但是它没有任何形状或者坐标信息。

以下为具有一个多边形、一个缺失值和一个空多边形的GeoSeries示例：

from shapely.geometry import Polygon
s = gpd.GeoSeries([Polygon([(0, 0), (1, 1), (0, 1)]), None, Polygon([])])
s

0    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....
1                                                 None
2                                        POLYGON EMPTY
dtype: geometry

在geopandas空间运算中，缺失的几何图形通常会传播。在结果中，这些缺失的几何图形也会缺失。另一方面，空的几何图形被视为几何图形。结果将取决于所进行的运算。如下所示：

s.area

0    0.5
1    NaN
2    0.0
dtype: float64

我们可以通过isna函数和is_empty属性判断是否为缺失值或者空值：

# 判断缺失值
s.isna()

0    False
1     True
2    False
dtype: bool

# 判断空值
www.gyballet.com_empty

0    False
1    False
2     True
dtype: bool

# 判断缺失或为空
www.gyballet.com_empty | s.isna()

0    False
1     True
2     True
dtype: bool

# 提取既不缺失也不为空的值
s[~(www.gyballet.com_empty | s.isna())]

0    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....
dtype: geometry

5 背景地图叠加

contextily是一个Python库，它提供了一种简单的方法将背景地图（通常是Web瓦片地图，如OpenStreetMap、Stamen Maps、Mapbox等）添加到地理空间数据可视化中。使用contextily库可以使地理空间数据可视化更加生动、直观，同时可以提供更多的地理信息。瓦片地图是一种基于网格的地图显示方式，将地图划分为多个小块，每个小块称为“瓦片”，每个瓦片都有自己的坐标和编号。这些瓦片可以按需加载，使用户能够快速地浏览地图，同时减少了加载时间和资源消耗。瓦片地图常用于在线地图应用程序，例如谷歌地图和百度地图。
contextily支持使用WGS84 (EPSG:4326)和Spheric Mercator (EPSG:3857)坐标系，在Web地图应用程序中，一般使用EPSG:3857（以米为单位）来显示瓦片地图，并使用EPSG:4326（以经纬度为单位）来标记瓦片地图上的位置。

contextily库的主要功能包括：

• 从Web地图提供商获取地图图层
• 将地图图层与地理空间数据集合并
• 使用Matplotlib或Bokeh绘制地图

本文主要介绍contextily简单使用，contextily具体使用可参考其官方文档：contextily-doc。contextily库中基于add_basemap函数在地图上添加背景地图。下面是该函数常用可用参数的介绍：

• ax: matplotlib axes对象，用于绘制地图
• crs: 输出地图的坐标系，默认为'EPSG:3857'
• source: 底图的来源，支持多种来源，如OpenStreetMap、Stamen Terrain、Stamen Toner等等，默认为OpenStreetMap
• zoom: 底图的缩放级别，默认为None，自动根据ax的extent和crs计算。zoom值越高，底图的缩放级别就越大，地图显示的范围也就越小，细节也会越来越清晰。
• url: 底图的url地址，默认为None，自动根据source和zoom计算。
• attribution: 底图的版权信息，默认为None
• alpha: 底图的透明度，默认为1.0
• *kwargs: 其他matplotlib.image()函数的可选参数，如cmap、vmin、vmax等等

⚠⚠source参数选择不同底图的来源，可能需要大量时间或者特定网络，如果失败多重试运行代码。

5.1 简单背景地图叠加

import geopandas as gpd
# 读取深圳市地图数据
data = www.gyballet.com_file("https://www.gyballet.com/areas_v3/bound/440300_full.json")
# 简单绘图
ax = data.plot(alpha=0.5, edgecolor="k")

# 确定数据所使用的坐标系
www.gyballet.com
# 将数据集所使用坐标系转为EPSG:3857
data_wm = www.gyballet.com_crs(epsg=3857)

import contextily as cx
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 10))
ax = data_wm.plot(ax = ax, alpha=0.5, edgecolor="k")
# 将自动下载瓦片地图
cx.add_basemap(ax)
# 保存地图
fig.savefig('save.jpg', pad_inches=0.5, bbox_inches='tight', dpi=300)

在上面的代码中，如果仅使用经纬度数据叠加瓦片地图，需要在add_basemap函数中设置crs参数，如下所示：

import contextily as cx
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 10))
ax = data.plot(ax = ax, alpha=0.5, edgecolor="k")
# 将自动下载瓦片地图
cx.add_basemap(ax, crs=www.gyballet.com)

可以通过调整zoom参数改变背景瓦片地图的细节程度，建议zoom值不要过大，下载速度太慢。此外可以通过设置attribution=""去除绘图水印。

ax = data_wm.plot(figsize=(10, 10), alpha=0.5, edgecolor="k")
cx.add_basemap(ax, zoom=12, attribution="")

5.2 定制化背景地图

通过设置add_basemap的source参数能够指定不同的数据源，以在地图上添加不同类型的底图。如下所示：

ax = data_wm.plot(figsize=(10, 10), alpha=0.5, edgecolor="k")
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite)
ax.set_axis_off()

当然也可以叠加多个背景地图，如下所示：

ax = data_wm.plot(figsize=(10, 10), alpha=0.5, edgecolor="k")
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite)
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels)

此外也可以叠加多个不同来源的背景图层。如下所示：

ax = data_wm.plot(figsize=(10, 10), alpha=0.5, edgecolor="k")
cx.add_basemap(ax, source=cx.providers.Stamen.Watercolor, zoom=12)
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels, zoom=10)

从上面的案例我们可以看到，contextily通过Provider预置提供商的名称来获取相应的Web瓦片地图。contextily所有预置的地图提供商可以通过以下cx.providers命令获取。可以尝试根据这些提供商定制瓦片地图格式。

# cx.providers

除了使用contextily预置的地图提供商，可以通过source设置给定瓦片地图地址来指定需要添加的底图。例如可以添加天地图，高德地图，腾讯地图的瓦片地图的地址。一些示例的瓦片地图地址可见：高德谷歌腾讯天地图地图瓦片url和在geopandas中叠加在线地图。
一般地图服务提供XYZ瓦片地图链接，其中的xyz代表了地图的坐标系。如下所示：

• x：表示在地图水平方向上的位置，从左到右递增，即经度值。
• y：表示在地图竖直方向上的位置，从上到下递增，即纬度值。
• z：表示地图的缩放级别，从0开始递增，数值越大，地图显示的范围越小，细节越丰富。

在瓦片地图中，地图被分成了许多小块，每个小块都有一个唯一的编号，也就是xyz坐标系。当我们使用地图服务时，通过改变xyz的值，就可以获取到不同位置、不同缩放级别下的地图瓦片，从而达到展示不同地图的目的。直接通过url设置瓦片地图示例如下：

fig, ax = plt.subplots(figsize=(10, 10))
ax = data_wm.plot(ax=ax, alpha=0.5, edgecolor="k")
cx.add_basemap(ax,
source='https://www.gyballet.com/appmaptile?style=6&x={x}&y={y}&z={z}',
zoom=12)
fig.savefig('save.jpg', pad_inches=0, bbox_inches='tight', dpi=300)

fig, ax = plt.subplots(figsize=(10, 10))
ax = data_wm.plot(ax=ax, alpha=0.5, edgecolor="k")
cx.add_basemap(ax,
source='https://www.gyballet.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
zoom=12)
fig.savefig('save.jpg', pad_inches=0, bbox_inches='tight', dpi=300)
ax.set_xlim(data_www.gyballet.com_bounds[0], data_www.gyballet.com_bounds[2])
ax.set_ylim(data_www.gyballet.com_bounds[1], data_www.gyballet.com_bounds[3])

(2559177.946084248, 2615308.057854809)

5.3 离线背景地图

在有些时候我们需要离线使用背景瓦片地图，contextly提供bounds2raster函数用于根据给定的空间范围和地图缩放级别，将在线地图服务中的栅格数据下载为本地文件。bounds2raster函数常用参数如下：

• w：float类型，表示空间范围的最小值。
• s：float类型，表示空间范围的最小值。
• e：float类型，表示空间范围的最大值。
• n：float类型，表示空间范围的最大值。
• path：str类型，表示下载的栅格数据文件的保存路径。
• zoom：int或者字符串类型，表示地图缩放级别。如果为字符串类型，可以设置为'auto'，表示自动确定最佳的缩放级别。
• source：str类型，表示地图服务的地址。
• ll：bool类型，表示w、s、e、n是否使用经纬度坐标系，默认为False。
• wait：int类型，表示两次下载之间的等待时间，单位为秒。默认为0。
• max_retries：int类型，表示下载失败后最大的重试次数。默认为2次。

bounds2raster函数返回RGB图像数组和瓦片图像边界框[minX，maxX，minY，maxY]，此外由于网络地图总是基于WGS84 Web Mercator（EPSG:3857）坐标系，因此bounds2raster函数返回和保存的图片都是基于EPSG:3857坐标系。

import geopandas as gpd
# 读取郑州市地图数据
data = www.gyballet.com_file("https://www.gyballet.com/areas_v3/bound/410100_full.json")
# 简单绘图
ax = data.plot(alpha=0.5, edgecolor="k")

# 叠加地图
ax = data.plot(alpha=0.5, edgecolor="k")
# crs告诉数据集用的坐标系统，这里www.gyballet.com为WGS 84(经纬度)
cx.add_basemap(ax,
crs=www.gyballet.com,
source=cx.providers.Stamen.Watercolor
)

提取待绘图区域的边界信息

west, south, east, north = bbox = www.gyballet.com_bounds
bbox

array([112.721178,  34.262109, 114.220962,  34.989506])

根据边界信息下载数据

import contextily as cx
import matplotlib.pyplot as plt
img, ext = cx.bounds2raster(west,
south,
east,
north,
"demo.tif",
source=cx.providers.Stamen.Watercolor,
ll=True
)

# 展示下载的数据
plt.axis('off')
plt.imshow(img)
# 边界范围
ext

(12523442.714243276, 12719121.506653327, 4030983.1236470537, 4187526.157575096)

有了背景地图，add_basemap中的source函数设置文件路径地址就可以离线叠加地图。

ax = data.plot(alpha=0.5, edgecolor="k")
# crs告诉数据集用的坐标系统，这里www.gyballet.com为WGS 84(经纬度)
cx.add_basemap(ax,
crs=www.gyballet.com,
source="demo.tif",
)

6 参考

• GeoPandas
• GeoPandas-doc
• GeoPandas Examples Gallery
• Python绘制数据地图1-GeoPandas入门指北
• Python绘制数据地图2-GeoPandas地图可视化
• matplotlib-scalebar
• contextily
• contextily-doc
• 高德谷歌腾讯天地图地图瓦片url
• 在geopandas中叠加在线地图