使用 Python 绘制城市间专利合作申请数量网络图(二)

继续上次的内容，今天我们使用 Python 的 matplotlib 来绘制城市间专利合作申请数量网络图。R 语言版本使用的是 ggraph 包，这里我们使用 matplotlib 来实现类似的网络图绘制效果。

也可以补充学习系列课程「ggplot2 数据可视化」的网络图绘制课时：https://rstata.duanshu.com/#/brief/course/9d8e0cc791644376979d78edc73eb18a

使用 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)

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

py_pkgs <- c(
  "numpy", "pandas", "geopandas", "shapely", "matplotlib", "scipy"
)
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", "shapely", "matplotlib")
pkgs[pkgs$package %in% key_pkgs, c("package", "version")]

虚拟环境管理常用命令

virtualenv_list()

# 删除虚拟环境（当不再需要时）
# virtualenv_remove(".venv")

# 升级某个包
# virtualenv_install(".venv", packages = "geopandas", ignore_installed = TRUE)

---

数据读取与预处理

加载 Python 包与字体

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
import matplotlib.patches as mpatches
import matplotlib.patheffects as pe
from shapely.geometry import LineString
import warnings
warnings.filterwarnings('ignore')

# 加载中文字体
font_path = "LXGWWenKai-Regular.ttf"
try:
    cnfont = fm.FontProperties(fname=font_path)
    fm.fontManager.addfont(font_path)
    font_name = fm.FontProperties(fname=font_path).get_name()
    plt.rcParams['font.family'] = font_name
    plt.rcParams['axes.unicode_minus'] = False
    print(f"字体加载成功: {font_name}")
except:
    cnfont = fm.FontProperties(family="SimHei")
    print("使用默认中文字体 SimHei")

配色方案与投影参数

# 配色方案
COL_PAL = [
    "#fed439", "#709ae1", "#8a9197", "#d2af81", "#fd7446",
    "#d5e4a2", "#197ec0", "#f05c3b", "#46732e", "#71d0f5"
]

# 线条样式
LINE_STYLES = {
    "九段线": {"color": "black", "linewidth": 0.5},
    "海岸线": {"color": "#0055AA", "linewidth": 0.2},
    "小地图框格": {"color": "black", "linewidth": 0.2},
    "省份": {"color": "#4d4d4d", "linewidth": 0.3},
}

# Albers 等面积投影（中国地图标准投影）
MY_CRS = ("+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")

读取数据

df = pd.read_stata("2020年城市间各类型专利合作数量统计.dta")
df = df[["城市1", "城市2", "合作申请专利数量"]].copy()
df.columns = ["from", "to", "value"]
# 过滤三沙市
df = df[(df["from"] != "三沙市") & (df["to"] != "三沙市")].reset_index(drop=True)
print(f"数据行数: {len(df)}")
df

读取地图数据

import re

# 省级线条
provline = gpd.read_file(
    "chinaprov2021mini/chinaprov2021mini_line.shp", encoding="utf-8"
)
mask = provline["class"].apply(
    lambda x: not re.search("_", str(x)) if pd.notna(x) else True
)
mask &= ~provline["class"].isin(["胡焕庸线", "秦岭-淮河线"])
provlinemap = provline[mask].copy()

# 城市面状底图
citymap = gpd.read_file(
    "chinacity2021mini/chinacity2021mini.shp", encoding="utf-8"
)
citymap = citymap[citymap["省代码"].notna()].copy()

# 完整版城市地图（计算质心用）
city = gpd.read_file("2021行政区划/市.shp", encoding="utf-8")
city = city.to_crs(MY_CRS)

print(f"省级线条: {len(provlinemap)} 条")
print(f"城市面: {len(citymap)} 个")

计算城市质心

# 使用 representative_point()（等价于 R 的 st_point_on_surface）
centroid_gdf = city[["省", "市", "geometry"]].copy()
centroid_gdf["geometry"] = city.geometry.representative_point()
centroid_df = centroid_gdf.copy()
centroid_df["X"] = centroid_gdf.geometry.x
centroid_df["Y"] = centroid_gdf.geometry.y
centroid_df = centroid_df[["省", "市", "X", "Y"]]
centroid_df = centroid_df[centroid_df["市"] != "三沙市"].reset_index(drop=True)
print(f"城市质心: {len(centroid_df)} 个")
centroid_df

网络数据构建

计算每个城市出发的总合作量

city_sum = df.groupby("from", as_index=False)["value"].sum()
city_sum.columns = ["from", "sum"]
city_sum = city_sum.sort_values("sum", ascending=False).reset_index(drop=True)
print("Top 10 城市:")
city_sum.head(10)

按 sum 排序并准备绘制数据

# 给 df 添加 sum 列并按 sum 排序
sum_map = dict(zip(city_sum["from"], city_sum["sum"]))
df["sum"] = df["from"].map(sum_map)
# 按升序排列：高 sum 的边绘制在上层（后绘制）
df = df.sort_values("sum").reset_index(drop=True)

# 按 sum 降序的 DataFrame（用于 Top N 选择）
sumdf = city_sum.copy()

# 建立城市到坐标的映射
coord_map = dict(zip(centroid_df["市"], zip(centroid_df["X"], centroid_df["Y"])))

# 简化地图数据
city_simplified = city.copy()
city_simplified = city_simplified.simplify(tolerance=2000, preserve_topology=True)
print("数据准备完成")

地图装饰函数

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

绘制简单网络图（直线版）

首先绘制最简单的直线网络图，每条边直接连接出发城市和目标城市的质心：

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

# 收集所有坐标确定范围（包含底图 bounds，确保小地图框格完整）
_map_bounds = citymap.total_bounds
all_x = [_map_bounds[0], _map_bounds[2]]
all_y = [_map_bounds[1], _map_bounds[3]]
for _, row in df.iterrows():
    if row["from"] in coord_map and row["to"] in coord_map:
        x0, y0 = coord_map[row["from"]]
        x1, y1 = coord_map[row["to"]]
        all_x.extend([x0, x1])
        all_y.extend([y0, y1])

xmin, xmax = min(all_x), max(all_x)
ymin, ymax = min(all_y), max(all_y)
margin_x = (xmax - xmin) * 0.05
margin_y = (ymax - ymin) * 0.05
xlim = (xmin - margin_x, xmax + margin_x)
ylim = (ymin - margin_y, ymax + margin_y)
ax.set_xlim(xlim)
ax.set_ylim(ylim)

# 底图
add_map_decorations(ax, citymap, provlinemap, cnfont, xlim, ylim)

# 绘制边
for _, row in df.iterrows():
    if row["from"] in coord_map and row["to"] in coord_map:
        x0, y0 = coord_map[row["from"]]
        x1, y1 = coord_map[row["to"]]
        ax.plot([x0, x1], [y0, y1], color="black", linewidth=0.05, zorder=3)

# 等比例
ax.set_aspect("equal")
_xlim = ax.get_xlim()
_ylim = ax.get_ylim()
data_w = _xlim[1] - _xlim[0]
data_h = _ylim[1] - _ylim[0]
fig.set_size_inches(10, 10 * data_h / data_w)
ax.set_xlim(_xlim)
ax.set_ylim(_ylim)

ax.set_title("geom_edge_bundle_force, width = 0.05",
             fontsize=14, fontweight="bold", pad=20, fontproperties=cnfont)
ax.set_axis_off()
fig.patch.set_facecolor("white")

plt.tight_layout()
fig.savefig("pic2020a.png", dpi=300, bbox_inches="tight",
            facecolor="white", edgecolor="none")
plt.close()
print("已保存: pic2020a.png")

绘制路径捆绑网络图

使用基于 Force-Directed Edge Bundling (FDEB) 的路径捆绑算法，通过 KDTree 空间索引加速邻近搜索，将相邻的边吸引到一起，使网络图更加清晰：

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

绘制最小捆绑网络图

使用较小的捆绑强度和迭代次数，产生更轻度的捆绑效果：

print("开始最小捆绑计算...")
bundled_minimal = edge_bundle_fdeb(df, centroid_df, strength=0.05, n_iter=30,
                                   subdivision=6, k_neighbors=15)

fig, ax = plt.subplots(figsize=(10, 8.5))
ax.set_xlim(xlim)
ax.set_ylim(ylim)

add_map_decorations(ax, citymap, provlinemap, cnfont, xlim, ylim)

for pts in bundled_minimal.values():
    line = LineString(pts)
    x, y = line.xy
    ax.plot(x, y, color="black", linewidth=0.05, zorder=3)

ax.set_aspect("equal")
_xlim = ax.get_xlim(); _ylim = ax.get_ylim()
dw = _xlim[1] - _xlim[0]; dh = _ylim[1] - _ylim[0]
fig.set_size_inches(10, 10 * dh / dw)
ax.set_xlim(_xlim); ax.set_ylim(_ylim)

ax.set_title("geom_edge_bundle_minimal, width = 0.05, max_distortion = 10",
             fontsize=14, fontweight="bold", pad=20, fontproperties=cnfont)
ax.set_axis_off()
fig.patch.set_facecolor("white")

plt.tight_layout()
fig.savefig("pic2020c.png", dpi=300, bbox_inches="tight",
            facecolor="white", edgecolor="none")
plt.close()
print("已保存: pic2020c.png")

绘制详细版网络图

下面我们复用图3的路径绑定结果（bundled_minimal），给图表添加更多细节——Top 100 城市节点、Top 10 城市标签、彩色边等：

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

这样我们就用 Python 绘制好了这幅地图，感兴趣的小伙伴也可以试试区县和省份的。

点击这里跳转到 RStata 短书平台获取附件：使用 Python 绘制城市间专利合作申请数量网络图(二)