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

前不久给大家分享了「1985～2024年省份间、城市间及区县间各类型专利合作数量面板数据」，在推文中我给大家展示了区县间、城市间和省份间的专利合作网络。今天我们以城市间专利申请合作网络为例讲解如何使用 Python 绘制这周空间网络图。

本次课程使用 Python 的 matplotlib + geopandas 复现了 R 语言中使用 ggplot2 基础图层绘制空间网络图的方法。

在学习本课程之前需要预先学习之前的课程：

使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）: https://rstata.duanshu.com/#/brief/course/379d19770956478ebc4d919910fca74c

使用 R 语言绘制城市间专利合作申请数量网络图（一）: https://rstata.duanshu.com/#/brief/course/xxxxx

本课程中使用的地图数据也都是来自该课程。

---

使用 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（仅首次）

# 设置 CRAN 镜像（knit 时 R 处于非交互模式，不会自动选择镜像）
options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

# 仅在尚未安装时才安装，避免每次 knit 都重装
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)
}

# 通过环境变量抢先锁定 Python（优先级最高，早于任何 {python} chunk）
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", "adjustText"
)
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 包已就绪，无需安装")
}

验证激活状态

# 验证当前绑定的 Python 路径（应指向 .venv 目录）
py_config()

查看已安装的包

# 列出虚拟环境中已安装的包
pkgs <- py_list_packages(".venv")
# 只显示我们关心的包
key_pkgs <- c("numpy", "pandas", "geopandas", "matplotlib", "shapely", "adjustText")
pkgs[pkgs$package %in% key_pkgs, c("package", "version")]

虚拟环境管理常用命令

# 查看所有已创建的虚拟环境
virtualenv_list()

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

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

---

数据读取与预处理

加载 Python 包

import pandas as pd
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
import matplotlib.patches as mpatches
from matplotlib.lines import Line2D
from matplotlib.colors import to_rgba
import os

# 加载中文字体
font_path = "LXGWWenKai-Regular.ttf"
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

设置坐标系

中国地图通常使用 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”

关键函数对照：

R	Python	说明
sf::st_transform(mycrs)	gdf.to_crs(MY_CRS)	坐标转换
sf::st_point_on_surface()	gdf.geometry.representative_point()	计算表面代表点
sf::st_coordinates()	gdf.geometry.x / .y	提取坐标
sf::st_simplify(dTolerance=2000)	gdf.geometry.simplify(tolerance=2000)	简化几何

读取专利合作数据

读取 2020 年城市间各类型专利合作数量统计数据：

df = pd.read_stata("2020年城市间各类型专利合作数量统计.dta")
df

数据清洗

重命名变量，并过滤掉三沙市（位置偏远且专利申请量极少）：

df = df.rename(columns={"城市1": "from", "城市2": "to",
                        "合作申请专利数量": "value"})
df = df[~df["from"].isin(["三沙市"]) & ~df["to"].isin(["三沙市"])]
df = df.reset_index(drop=True)

关键函数对照：

R (dplyr)	Python (pandas)	说明
rename(from = 城市1, to = 城市2)	rename(columns={"城市1": "from", ...})	重命名列
filter(!from %in% "三沙市")	df[~df["from"].isin(["三沙市"])]	条件筛选
ungroup()	无需操作	pandas 无 group 状态

统计每个城市的总合作量

# 每个城市的总合作量（等价于 R 的 group_by + summarise）
dfsum = df.groupby("from")["value"].sum().reset_index()
dfsum.columns = ["from", "value"]
dfsum

# 每条记录所在城市的总合作量（等价于 R 的 group_by + mutate(sum = sum(value))）
city_total = df.groupby("from")["value"].sum().reset_index()
city_total.columns = ["from", "sum"]
df = df.merge(city_total, on="from", how="left")

---

地图数据准备

读取地图数据

# 省界线（九段线、海岸线等）
provlinemap = gpd.read_file(
    "chinaprov2021mini/chinaprov2021mini_line.shp", encoding="utf-8"
)
provlinemap = provlinemap[
    ~provlinemap["class"].str.contains("_", na=False) &
    ~provlinemap["class"].isin(["胡焕庸线", "秦岭-淮河线"])
]
provlinemap = provlinemap[["class", "geometry"]]

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

关键函数对照：

R (sf)	Python (geopandas)	说明
read_sf("xxx.shp")	gpd.read_file("xxx.shp")	读取矢量
filter(!str_detect(class, "_"))	~df["class"].str.contains("_")	正则筛选
select(class)	df[["class", "geometry"]]	选择列
filter(!is.na(省代码))	df[df["省代码"].notna()]	去除空值

计算城市质心

读取未经编辑的城市矢量数据，计算各城市的质心坐标（等价于 R 的 st_point_on_surface）：

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

匹配合作线起点和终点坐标

df = df.merge(city_centroiddf, left_on="from", right_on="市", how="left")
df = df.rename(columns={"X": "fromlon", "Y": "fromlat", "省": "prov"})
df = df.drop(columns=["市"])

df = df.merge(city_centroiddf, left_on="to", right_on="市", how="left")
df = df.rename(columns={"X": "tolon", "Y": "tolat"})
df = df.drop(columns=["市", "省"], errors="ignore")

df

简化地图数据

city_simplified = city.geometry.simplify(tolerance=2000, preserve_topology=True)
city = city.set_geometry(city_simplified)

---

准备散点和标签数据

由于城市数量众多，需要筛选出最重要的城市进行显示：

• 前 100 个城市绘制散点（按总合作量排序）
• 前 10 个城市添加文本标签（避免标签过于密集）

# 前100个城市用于散点（等价于 R 的 slice(1:100)）
first100df = (
    df.groupby("from")
    .agg(sum=("value", "sum"),
         fromlon=("fromlon", "first"),
         fromlat=("fromlat", "first"))
    .reset_index()
    .sort_values("sum", ascending=False)
    .head(100)
    .reset_index(drop=True)
)

# 前10个城市用于标签（等价于 R 的 slice(1:10)）
first10df = first100df.head(10).copy()
first10df

---

绘制网络图

配色方案

地图配色可以使用这个网站：https://tidyfriday.cn/colors

# 线要素颜色和线宽
line_colors = {
    "九段线": "black",
    "海岸线": "#0055AA",
    "小地图框格": "black",
    "省份": "#b3b3b3"  # R 的 gray30
}
line_widths = {
    "九段线": 0.5,
    "海岸线": 0.2,
    "小地图框格": 0.2,
    "省份": 0.3
}

# 前10个高亮城市的颜色
top10_colors = ["#fed439", "#709ae1", "#8a9197", "#d2af81", "#fd7446",
                "#d5e4a2", "#197ec0", "#f05c3b", "#46732e", "#71d0f5"]

绘图核心代码

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

R vs Python 关键函数对照

本次课程涉及的 R → Python 核心函数映射：

R (ggplot2/sf)	Python (matplotlib/geopandas)	说明
geom_sf(fill=NA, color="gray10")	gdf.plot(ax=ax, facecolor="none", edgecolor="#1a1a1a")	绘制面边界
geom_curve(curvature=0.3)	ax.annotate(..., connectionstyle="arc3,rad=0.3")	弧线连接
scale_linewidth_continuous(range=c(0.1, 0.5))	0.1 + 0.4 * (x - min) / (max - min)	线宽映射
scale_size_continuous(range=c(0.01, 7))	0.01 + 6.99 * (x - min) / (max - min)	散点大小映射
scale_alpha_manual(values=c(rep(1,10), rep(0.2, n-10)))	alpha = 1.0 if city in top10 else 0.2	透明度控制
ggrepel::geom_text_repel()	adjustText.adjust_text()	避免重叠文本标签
annotation_scale()	手绘比例尺	比例尺
annotation_north_arrow()	ax.annotate(..., arrowprops=dict(arrowstyle="->"))	指北针
guide_legend(direction="horizontal")	ax.legend(ncol=2)	图例布局
st_point_on_surface()	gdf.geometry.representative_point()	表面代表点
st_coordinates()	gdf.geometry.x / .y	提取坐标

---

省份间和区县间网络图

附件中也提供了省份和区县版本的绘图代码（main2.py 和 main3.py），感兴趣的小伙伴也可以拿来参考。主要的区别在于：

1. 数据源不同：省份版使用 2020年省份间各类型专利合作数量统计.dta，区县版使用 2020年区县间各类型专利合作数量统计.dta
2. ID 字段不同：省份版用 省1/省2 匹配省份名称，区县版用 县代码1/县代码2 匹配区县代码
3. 地图底图不同：省份版使用 chinaprov2021mini，区县版使用 chinacounty2021mini
4. 线宽范围不同：省份版使用 range=c(0.2, 1.5)，区县版使用 range=c(0.1, 0.5)
5. 透明度不同：区县版非 top10 城市使用 alpha=0.05（因为区县数量更多，需要更低的透明度）
6. 标签内容不同：区县版标签使用 市+县 格式，而非单纯的区县名称

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