使用 Python 绘制带邻国地区的中国地图

在之前「使用 ggplot2 绘制地图——以中国地图为例」课程的基础上，我们今天再来学习下如何使用 Python 的 matplotlib 把邻国的区域也添加到中国地图上：

小地图版本	长版

使用 reticulate 创建与管理 Python 虚拟环境

在 R 中通过 reticulate 包来调用 Python，最好的实践是为项目创建一个专属的 Python 虚拟环境，将所需依赖隔离到独立空间，避免与系统 Python（如 Anaconda）发生版本冲突。

重要说明（避免”已初始化”报错）：reticulate 在 R 会话中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定，之后再调用 use_virtualenv() 会报错：

ERROR: The requested version of Python cannot be used, as another version has already been initialized.

因此，虚拟环境的激活必须在所有 {python} 代码块之前完成。本文档的解决方案是在 setup chunk 中通过 Sys.setenv(RETICULATE_PYTHON = ...) 提前锁定 Python 路径，这是 reticulate 选取 Python 的最高优先级入口。

安装 reticulate（仅首次）

options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

if (!requireNamespace("reticulate", quietly = TRUE)) {
  install.packages("reticulate")
  message("reticulate 安装完成！")
} else {
  message("reticulate 已安装，版本：", packageVersion("reticulate"))
}

虚拟环境初始化原理（已在 setup chunk 中完成）

本文档的 setup chunk（隐藏运行）包含如下逻辑：

library(reticulate)

.venv_name   <- ".venv"
.venv_python <- virtualenv_python(.venv_name)

if (!file.exists(.venv_python)) {
  virtualenv_create(.venv_name)
  .venv_python <- virtualenv_python(.venv_name)
}

Sys.setenv(RETICULATE_PYTHON = .venv_python)
use_virtualenv(.venv_name, required = TRUE)

这样做的关键在于：knitr 在处理第一个 {python} chunk 时，reticulate 已经通过 RETICULATE_PYTHON 环境变量知道要使用 .venv，不会再去碰 Anaconda。

在虚拟环境中安装 Python 包（仅首次）

py_pkgs <- c(
  "numpy", "pandas", "geopandas", "matplotlib",
  "shapely", "openpyxl"
)
installed <- py_list_packages(".venv")$package
need_install <- setdiff(py_pkgs, installed)

if (length(need_install) > 0) {
  virtualenv_install(".venv", packages = need_install)
  message("已安装缺失的包：", paste(need_install, collapse = ", "))
} else {
  message("所有 Python 包已就绪，无需安装")
}

验证激活状态

py_config()

查看已安装的包

pkgs <- py_list_packages(".venv")
key_pkgs <- c("numpy", "pandas", "geopandas", "matplotlib", "shapely", "openpyxl")
pkgs[pkgs$package %in% key_pkgs, c("package", "version")]

虚拟环境管理常用命令

virtualenv_list()
# virtualenv_remove(".venv")  # 删除虚拟环境
# virtualenv_install(".venv", packages = "geopandas", ignore_installed = TRUE)  # 升级某个包

---

小地图版本

1. 在地图上添加散点

首先加载所需 Python 包，设置字体和坐标系：

import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.font_manager as fm
from matplotlib.patches import Rectangle
from shapely.geometry import box, Point
import warnings
warnings.filterwarnings("ignore")

# 字体设置
font_path = "LXGWWenKai-Regular.ttf"
font_prop = fm.FontProperties(fname=font_path)
fm.fontManager.addfont(font_path)
plt.rcParams["font.family"] = font_prop.get_name()
plt.rcParams["axes.unicode_minus"] = False

# 中国地图投影坐标系
MYCRS = (
    "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 "
    "+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
)

然后读取地图数据和散点数据：

# 绘图区域（预设好的 bbox）
plot_bbox = box(-2725586, 1800655, 2982768, 6000000)
plot_bbox_gdf = gpd.GeoDataFrame(geometry=[plot_bbox], crs=MYCRS)

# 读取小地图版本的中国城市地图
citymap = gpd.read_file("chinacity2019mini/chinacity2019mini.shp")
citymap = citymap[citymap["省代码"].notna()].copy()
citymap = citymap.to_crs(MYCRS)

# 线条（九段线、海岸线、小地图框格）
citylinemap_all = gpd.read_file("chinacity2019mini/chinacity2019mini_line.shp")
citylinemap_all = citylinemap_all.to_crs(MYCRS)
citylinemap = citylinemap_all[
    citylinemap_all["class"].isin(["九段线", "海岸线", "小地图框格"])
].copy()

在之前课程的基础上，我们还需要读取一份邻国地图数据：

# 邻国地图
china_neighboring = gpd.read_file("china_neighboring/china_neighboring.shp")
china_neighboring = china_neighboring.to_crs(MYCRS)

读取企业散点数据，并将经纬度坐标转换为投影坐标：

# 读取瞪羚、独角兽、创新型企业经纬度数据
pointdf = pd.read_excel("瞪羚、独角兽、创新型企业经纬度数据.xlsx")
pointdf = pointdf[pointdf["经度"].notna()].copy()
print(f"数据量：{len(pointdf)} 条")

# 转换成 GeoDataFrame（WGS84 → 投影坐标）
pointdfsf = gpd.GeoDataFrame(
    pointdf,
    geometry=gpd.points_from_xy(pointdf["经度"], pointdf["纬度"]),
    crs="EPSG:4326"
)
pointdfsf = pointdfsf.to_crs(MYCRS)

接下来处理小地图框格内的散点——先提取小地图范围内的点，然后将其缩放并平移到小地图框格的位置：

# 小地图范围
small_bbox = box(120000, 320000, 1766004.1, 2557786.0)

# 提取小地图范围内的点
pointdfsf_small = pointdfsf[pointdfsf.geometry.within(small_bbox)].copy()
print(f"小地图内散点数：{len(pointdfsf_small)}")

# 缩放 0.5 倍 + 平移（等价于 R 中 geometry * 0.5 + c(2100000, 1665139)）
new_geoms = []
for geom in pointdfsf_small.geometry:
    new_geoms.append(Point(geom.x * 0.5 + 2100000, geom.y * 0.5 + 1665139))

pointdfsf_small = pointdfsf_small.copy()
pointdfsf_small = pointdfsf_small.set_geometry(new_geoms)
pointdfsf_small = pointdfsf_small.set_crs(MYCRS)

# 合并主地图散点 + 小地图散点
pointdfsfall = pd.concat([pointdfsf, pointdfsf_small], ignore_index=True)
pointdfsfall = gpd.GeoDataFrame(pointdfsfall, geometry="geometry", crs=MYCRS)

设置省份颜色表（对应 R 中 ggsci::default_igv）和线条样式：

# 省份颜色表（32 色，对应 ggsci::default_igv）
IGV_COLORS = [
    "#5050FF", "#CE3D32", "#749B58", "#F0E685", "#466983", "#B7DEE8",
    "#70C4BF", "#56B4E9", "#E69F00", "#009E73", "#F0E442", "#0072B2",
    "#D55E00", "#CC79A7", "#999999", "#FF7F00", "#A6CEE3", "#1F78B4",
    "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", "#FDBF6F", "#FF7F00",
    "#CAB2D6", "#6A3D9A", "#FFFF99", "#B15928", "#FBB4AE", "#B3CDE3",
    "#CCEBC5", "#DECBE4",
]
province_list = sorted(pointdfsfall["省"].dropna().unique())
color_map = {prov: IGV_COLORS[i % len(IGV_COLORS)] for i, prov in enumerate(province_list)}
pointdfsfall["color"] = pointdfsfall["省"].map(color_map)

# 线条颜色与线宽
LINE_COLOR = {"九段线": "#A29AC4", "海岸线": "#0055AA", "小地图框格": "black"}
LINE_WIDTH = {"九段线": 0.6, "海岸线": 0.3, "小地图框格": 0.3}

然后就可以绘制散点地图了：

# 定义比例尺和指北针辅助函数
def add_scale_bar(ax, ref_gdf, location="bl"):
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    ax_w = xlim[1] - xlim[0]
    ax_h = ylim[1] - ylim[0]
    # 固定 1000km
    bar_km = 1000
    bar_m = bar_km * 1000
    bar_h = ax_h * 0.016
    bar_x = xlim[0] + ax_w * 0.03
    bar_y = ylim[0] + ax_h * 0.05
    bar_w = bar_m
    segments = 4
    seg_w = bar_w / segments
    for i in range(segments):
        fc = "black" if i % 2 == 0 else "white"
        rect = Rectangle(
            (bar_x + i * seg_w, bar_y), seg_w, bar_h,
            facecolor=fc, edgecolor="black", linewidth=0.6,
            clip_on=False, zorder=10
        )
        ax.add_patch(rect)
    ax.text(bar_x, bar_y - bar_h * 0.5, "0",
            fontsize=9, ha="center", va="top", fontproperties=font_prop, zorder=11)
    ax.text(bar_x + bar_w, bar_y - bar_h * 0.5, f"{bar_km} km",
            fontsize=9, ha="center", va="top", fontproperties=font_prop, zorder=11)

def add_north_arrow(ax):
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    ax_w = xlim[1] - xlim[0]
    ax_h = ylim[1] - ylim[0]
    x_center = xlim[1] - ax_w * 0.05
    y_base = ylim[1] - ax_h * 0.08
    arrow_len = ax_h * 0.06
    ax.annotate(
        "", xy=(x_center, y_base + arrow_len), xytext=(x_center, y_base),
        arrowprops=dict(arrowstyle="-|>", color="black", lw=2), zorder=11
    )
    ax.text(x_center, y_base + arrow_len + ax_h * 0.015, "N",
            fontsize=12, ha="center", va="bottom",
            fontweight="bold", fontproperties=font_prop, zorder=11)

此处代码需下载讲义材料查看～

2. 填充地图

下面我们再演示下填充地图的绘制。

首先统计每个城市的公司数量，并与地图数据合并：

citydf = pointdf.groupby(["市", "市代码"]).size().reset_index(name="n")
citymap2 = citymap.merge(citydf, on=["市", "市代码"], how="left")
citymap2["n"] = citymap2["n"].fillna(0).astype(int)
print(citymap2["n"].describe())

然后绘制连续填充地图：

n = n + 1 是为了对数化之后 0 变成 -Inf，在图上会反映为缺失值。

此处代码需下载讲义材料查看～

也可以使用分段填色：

# 计算分位数断点
cutlist_raw = np.quantile(citymap2["n"].values, np.arange(1, 11) / 10)
cutlist = sorted(set(cutlist_raw.astype(int).tolist()))

labels_c = ["<= 1", "1～3", "3～9", "9～18", "19～39", "39～92", ">= 92"]

citymap3 = citymap2.copy()
citymap3["group"] = pd.cut(
    citymap3["n"], bins=cutlist,
    labels=labels_c[: len(cutlist) - 1],
    include_lowest=True
)
print(citymap3["group"].value_counts().sort_index())

# 离散色阶（使用 inferno_r 的 0~0.95 范围）
n_cats = len(labels_c)
acton_colors = [CMAP_ACTON(i / (n_cats - 1) * 0.95) for i in range(n_cats)]
cat_color_map = dict(zip(labels_c, acton_colors))

fig, ax = plt.subplots(figsize=(11.5, 8.5))

plot_bbox_gdf.plot(ax=ax, color="#BBD1EB", zorder=1)

# 逐组绘制省份
for label in labels_c:
    subset = citymap3[citymap3["group"] == label]
    if len(subset) > 0:
        subset.plot(
            ax=ax, facecolor=cat_color_map[label],
            edgecolor="gray", linewidth=0.01, zorder=2
        )

# 未分组（NaN）用白色
citymap3[citymap3["group"].isna()].plot(
    ax=ax, facecolor="white", edgecolor="gray", linewidth=0.01, zorder=2
)

china_neighboring.plot(ax=ax, facecolor="#EDEDED", edgecolor="none", zorder=3)

for _, row in china_neighboring.iterrows():
    centroid = row.geometry.centroid
    ax.text(
        centroid.x, centroid.y, row["country_cn"],
        fontsize=8, color="gray", ha="center", va="center",
        fontproperties=font_prop, zorder=6
    )

for cls in ["九段线", "海岸线", "小地图框格"]:
    subset = citylinemap[citylinemap["class"] == cls]
    if len(subset) > 0:
        subset.plot(ax=ax, color=LINE_COLOR[cls], linewidth=LINE_WIDTH[cls], zorder=4)

ax.scatter(
    pointdfsfall.geometry.x, pointdfsfall.geometry.y,
    c=pointdfsfall["color"], s=0.3, linewidths=0, zorder=5
)

ax.set_aspect("equal")
ax.set_xlim(plot_bbox_gdf.total_bounds[0], plot_bbox_gdf.total_bounds[2])
ax.set_ylim(plot_bbox_gdf.total_bounds[1], plot_bbox_gdf.total_bounds[3])

xlim = ax.get_xlim()
ylim = ax.get_ylim()
ax_w = xlim[1] - xlim[0]
ax_h = ylim[1] - ylim[0]

# 手动绘制分类图例（左下角，两列排列，从上到下、从左到右）
legend_x = xlim[0] + ax_w * 0.03
legend_y = ylim[0] + ax_h * 0.10
sq_size = ax_h * 0.032
gap = ax_h * 0.040
n_cols = 2
col_offset = ax_w * 0.15
n_rows_per_col = (len(labels_c) + n_cols - 1) // n_cols

此处代码需下载讲义材料查看～

长版地图

长版地图的绘制更简单一些，这里演示省级分段填充地图的绘制。

首先读取长版省级地图、线条和邻国数据：

# 长版绘图区域
long_bbox = box(-2925762.0, 377031.1, 2507277.8, 6221888.6)
long_bbox_gdf = gpd.GeoDataFrame(geometry=[long_bbox], crs=MYCRS)

# 读取长版省级地图
provmap = gpd.read_file("chinaprov2021long/chinaprov2021long.shp")
provmap = provmap[provmap["省代码"].notna()].copy()
provmap = provmap.to_crs(MYCRS)

# 线条（长版无小地图框格）
provlinemap_all = gpd.read_file("chinaprov2021long/chinaprov2021long_line.shp")
provlinemap_all = provlinemap_all.to_crs(MYCRS)
provlinemap = provlinemap_all[
    provlinemap_all["class"].isin(["九段线", "海岸线"])
].copy()

# 邻国地图（长版）
china_neighboring_long = gpd.read_file("china_neighboring_long/china_neighboring_long.shp")
china_neighboring_long = china_neighboring_long.to_crs(MYCRS)

读取散点数据并汇总统计每个省份的企业数量：

# 读取散点数据
pointdf = pd.read_excel("瞪羚、独角兽、创新型企业经纬度数据.xlsx")
pointdf = pointdf[pointdf["经度"].notna()].copy()

pointdfsf = gpd.GeoDataFrame(
    pointdf,
    geometry=gpd.points_from_xy(pointdf["经度"], pointdf["纬度"]),
    crs="EPSG:4326"
)
pointdfsf = pointdfsf.to_crs(MYCRS)

# 省份颜色
province_list_long = sorted(pointdfsf["省"].dropna().unique())
color_map_long = {prov: IGV_COLORS[i % len(IGV_COLORS)] for i, prov in enumerate(province_list_long)}
pointdfsf["color"] = pointdfsf["省"].map(color_map_long)

# 汇总每个省份的企业数量
provdf = pointdfsf.drop(columns="geometry").groupby(["省", "省代码"]).size().reset_index(name="n")
provdf2 = provmap.merge(provdf, on=["省", "省代码"], how="left")
provdf2["n"] = provdf2["n"].fillna(0).astype(int)

按分位数分组：

cutlist_raw = np.quantile(provdf2["n"].values, np.arange(1, 9) / 8)
cutlist = [0] + sorted(set(cutlist_raw.astype(int).tolist()))

labels_p = ["<= 8", "8～28", "28～58", "58～159",
            "159～513", "513～870", "870～2459", "> 2450"]

provdf3 = provdf2.copy()
provdf3["group"] = pd.cut(
    provdf3["n"], bins=cutlist,
    labels=labels_p[: len(cutlist) - 1],
    include_lowest=True
)
print(provdf3["group"].value_counts().sort_index())

然后绘制长版省级填充地图：

# 长版线条样式
LINE_COLOR_LONG = {"九段线": "#A29AC4", "海岸线": "#0055AA"}
LINE_WIDTH_LONG = {"九段线": 0.6, "海岸线": 0.3}

# 离散色阶
n_cats = len(labels_p)
acton_colors_long = [CMAP_ACTON(i / (n_cats - 1) * 0.95) for i in range(n_cats)]
cat_color_map_long = dict(zip(labels_p, acton_colors_long))

fig, ax = plt.subplots(figsize=(8, 9))

long_bbox_gdf.plot(ax=ax, color="#BBD1EB", zorder=1)

# 逐组绘制省份
for label in labels_p:
    subset = provdf3[provdf3["group"] == label]
    if len(subset) > 0:
        subset.plot(
            ax=ax, facecolor=cat_color_map_long[label],
            edgecolor="gray", linewidth=0.01, zorder=2
        )

provdf3[provdf3["group"].isna()].plot(
    ax=ax, facecolor="white", edgecolor="gray", linewidth=0.01, zorder=2
)

china_neighboring_long.plot(ax=ax, facecolor="#EDEDED", edgecolor="none", zorder=3)

for _, row in china_neighboring_long.iterrows():
    centroid = row.geometry.centroid
    ax.text(
        centroid.x, centroid.y, row["country_cn"],
        fontsize=7, color="gray", ha="center", va="center",
        fontproperties=font_prop, zorder=6
    )

for cls in ["九段线", "海岸线"]:
    subset = provlinemap[provlinemap["class"] == cls]
    if len(subset) > 0:
        subset.plot(ax=ax, color=LINE_COLOR_LONG[cls], linewidth=LINE_WIDTH_LONG[cls], zorder=4)

ax.scatter(
    pointdfsf.geometry.x, pointdfsf.geometry.y,
    c=pointdfsf["color"], s=0.3, linewidths=0, zorder=5
)

ax.set_aspect("equal")
ax.set_xlim(long_bbox_gdf.total_bounds[0], long_bbox_gdf.total_bounds[2])
ax.set_ylim(long_bbox_gdf.total_bounds[1], long_bbox_gdf.total_bounds[3])

xlim = ax.get_xlim()
ylim = ax.get_ylim()
ax_w = xlim[1] - xlim[0]
ax_h = ylim[1] - ylim[0]

# 分类图例（左侧，下移到 0.12 位置，从上到下排列）
此处代码需下载讲义材料查看～

关键函数对照表

R（ggplot2 + sf）	Python（matplotlib + geopandas）	说明
st_bbox() %>% st_as_sfc()	shapely.geometry.box()	创建矩形区域
read_sf()	gpd.read_file()	读取 Shapefile
st_as_sf(coords, crs)	gpd.points_from_xy() + GeoDataFrame(crs)	创建点要素
st_transform(crs)	gdf.to_crs(crs)	坐标系转换
st_intersection()	gdf.geometry.within()	空间裁剪（点在面内）
bind_rows()	pd.concat()	合并数据框
geom_sf()	gdf.plot()	绘制地图图层
annotation_scale()	手绘黑白矩形比例尺	比例尺
annotation_north_arrow()	ax.annotate(arrowprops)	指北针
scale_fill_scico()	plt.cm.get_cmap() + mcolors.LogNorm()	连续色阶
scale_fill_scico_d()	pd.cut() + 逐组绘制	分段色阶
ggsave()	plt.savefig()	保存图片

点击这里跳转到 RStata 短书平台获取附件：使用 Python 绘制带邻国地区的中国地图