使用 Python 绘制城市间专利合作申请数量网络图(二)
- 2026-07-01 04:32:58
使用 Python 绘制城市间专利合作申请数量网络图(二)
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")) } 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) 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)
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) elseTrue ) 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) # 给 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 _, rowin df.iterrows(): if row["from"] in coord_map androw["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 _, rowin df.iterrows(): if row["from"] in coord_map androw["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") 
此处代码需下载讲义材料查看～ 
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") 
此处代码需下载讲义材料查看～ 
由于借助 AI 工具学习编程已经变得非常容易了，因此之后的课程就不再默认进行视频讲解了，如果特别需要视频讲解也可以联系李老师预约讲解～讲义材料学习过程中遇到的问题也可以及时与李老师联系。
继续上次的内容，今天我们使用 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}代码块之前完成。本文档的解决方案是在setupchunk 中通过Sys.setenv(RETICULATE_PYTHON = ...)提前锁定 Python 路径，这是 reticulate 选取 Python 的最高优先级入口。
安装 reticulate（仅首次）
虚拟环境初始化原理（已在 setup chunk 中完成）
本文档的 setup chunk（隐藏运行）包含如下逻辑：
在虚拟环境中安装 Python 包（仅首次）
验证激活状态
查看已安装的包
虚拟环境管理常用命令
数据读取与预处理
加载 Python 包与字体
配色方案与投影参数
读取数据
读取地图数据
计算城市质心
网络数据构建
计算每个城市出发的总合作量
按 sum 排序并准备绘制数据
地图装饰函数
绘制简单网络图（直线版）
首先绘制最简单的直线网络图，每条边直接连接出发城市和目标城市的质心：
绘制路径捆绑网络图
使用基于 Force-Directed Edge Bundling (FDEB) 的路径捆绑算法，通过 KDTree 空间索引加速邻近搜索，将相邻的边吸引到一起，使网络图更加清晰：
绘制最小捆绑网络图
使用较小的捆绑强度和迭代次数，产生更轻度的捆绑效果：
绘制详细版网络图
下面我们复用图3的路径绑定结果（bundled_minimal），给图表添加更多细节——Top 100 城市节点、Top 10 城市标签、彩色边等：
这样我们就用 Python 绘制好了这幅地图，感兴趣的小伙伴也可以试试区县和省份的。
如何参加课程？
是不是感觉很硬核！欢迎报名 RStata 培训班获取全部课程和以会员价获取数据资料（10元/份）详情可阅读这篇推文：数据处理、图表绘制、效率分析与计量经济学如何学习～
详情可点击阅读原文进入 RStata 学院了解（从首页的会员卡专区即可查看和购买会员卡）。
更多关于 RStata 培训班的信息可添加微信号 r_stata2 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/36e42f1c75fa482e80a078adf5f915e4
本文来自网友投稿或网络内容，如有侵犯您的权益请联系我们删除，联系邮箱：wyl860211@qq.com 。
