由于借助 AI 工具学习编程已经变得非常容易了，因此之后的课程就不再默认进行视频讲解了，如果特别需要视频讲解也可以联系李老师预约讲解～讲义材料学习过程中遇到的问题也可以及时与李老师联系。
购买 RStata 名师讲堂会员即可参加该课程啦（之前的和未来的都可以参加）！
价格：2800/年 或者 4800/长期
购买会员可以从这里下单：https://rstata.duanshu.com/#/card/list/
名师讲堂会员权益：
- 参加平台上的其他 R 语言和 Stata 的课程；
- 以会员折扣价购买我们分享的数据资料（10 元/份）；
* 如果发票可添加小编微信 r_stata2 （RStata 李老师）开具。如需数据资料，购买后可添加小编微信免费领取数据折扣卡。
更多关于 RStata 会员的更多信息可添加微信号 r_stata2 咨询：
课程主页（点击文末的阅读原文即可跳转）：https://rstata.duanshu.com/#/brief/course/eb9f29edaadb45728360ba8f6aa0c399
今天给大家介绍如何使用 Python 基于 LandScan 全球人口栅格数据，按照商玉萍（2022）论文中的方法，同时测算城市多中心结构的 5 大指标，包括中心数量、帕累托指数、多中心指数（含距离）以及去中心化指标。
一、指标来源与计算原理
1.1 数据来源
本方法使用的人口数据为美国能源部橡树岭国家实验室提供的 LandScan 全球人口密度栅格数据，空间分辨率约为 1 km × 1 km（坐标参考系转换后约 820 m）。
1.2 文献来源
5 大指标均来自以下论文：
商玉萍. 中国城市多中心空间战略的创新绩效研究——基于集聚经济与舒适度的视角 [J]. 经济学（季刊）, 2022.
该论文从集聚经济与舒适度双重视角，考察城市多中心空间战略对创新绩效的影响，所使用的多中心测量指标体系是目前文献中最为系统的之一。
1.3 五大指标定义
| | |
|---|
center | | |
pareto | | 各中心人口规模的秩-规模幂律回归系数（绝对值），越小表示多中心越均衡 |
poly | | |
sub3 | | |
sub5 | | |
1.4 各指标计算方法
1.5 计算流程总览
整个计算分为以下步骤：
- 坐标系转换：将 LandScan 栅格数据投影到等面积坐标系（Albers 投影），确保距离计算准确；
- 裁剪与掩膜：按城市行政边界裁剪栅格，获取该城市范围内的人口格点；
- 构建点邻接：使用 1000 m 距离阈值构建邻接关系，生成空间权重矩阵；
- 局部 Moran's I 检验：使用 R
spdep::localmoran(conditional=TRUE) 的条件正态近似法（Sokal 1998）计算每个格点的局部莫兰指数和 p 值，识别统计显著的 HH 聚集区格点； - 聚类成中心：将 HH 格点按空间邻接分组，形成连续的"高密度区块"；
- 筛选有效中心：要求每个区块格点数 ≥ 3、总人口 ≥ 100,000；
- 计算 5 大指标：center、pareto、poly、sub3、sub5。
1.6 参数设定说明
"sig_level": 0.05, # 局部莫兰指数显著性阈值"min_cells": 3, # 中心最小格点数（过滤噪点）"min_pop": 100000, # 中心最小人口（单位：人），论文基准为 10 万"dist_nb": 1000# 邻接距离阈值（米）for k, v in ANALYSIS_PARAMS.items():关于 dist_nb = 1000 的设定：LandScan 数据经 Albers 投影后分辨率约为 820 m。相邻格点（上下左右）的距离约为 820 m，对角线方向约为 820 × √2 ≈ 1159 m。将阈值设为 1000 m，可精确选取边邻接格点，不会错误地把对角线邻居算进来。
二、使用 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 的最高优先级入口。
2.1 安装 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 安装完成！") message("reticulate 已安装，版本：", packageVersion("reticulate"))2.2 虚拟环境初始化原理（已在 setup chunk 中完成）
本文档的 setup chunk（隐藏运行）包含如下逻辑：
.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。
2.3 在虚拟环境中安装 Python 包（仅首次）
"numpy", "pandas", "geopandas", "rasterio","scipy", "scikit-learn", "libpysal", "pyreadstat"installed <- py_list_packages(".venv")$packageneed_install <- setdiff(py_pkgs, installed)if (length(need_install) > 0) { virtualenv_install(".venv", packages = need_install) message("已安装缺失的包：", paste(need_install, collapse = ", ")) message("所有 Python 包已就绪，无需安装")2.4 验证激活状态
# 验证当前绑定的 Python 路径（应指向 .venv 目录）2.5 查看已安装的包
pkgs <- py_list_packages(".venv")key_pkgs <- c("numpy", "pandas", "geopandas", "rasterio", "libpysal", "scipy", "scikit-learn", "pyreadstat")pkgs[pkgs$package %in% key_pkgs, c("package", "version")]2.6 虚拟环境管理常用命令
# virtualenv_remove(".venv")# virtualenv_install(".venv", packages = "geopandas", ignore_installed = TRUE)
三、全局路径与参数配置
3.1 导入依赖包与参数设定
from rasterio.mask import maskfrom scipy.stats import normfrom sklearn.linear_model import LinearRegressionfrom libpysal.weights import DistanceBand, lag_spatialMYCRS = "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs""pop_tif": "pop-tif2", # 已重投影到 mycrs 的 tif 文件夹"city_shp": "2021行政区划/市.shp"# 2021 年市级行政区划数据预处理提示：如果你的 LandScan tif 文件还是 WGS84 坐标系，需要先转换到 MYCRS：
from rasterio.warp import calculate_default_transform, reproject, Resamplingdefreproject_raster(input_path, output_path, dst_crs):with rasterio.open(input_path) as src: transform, width, height = calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.boundswith rasterio.open(output_path, 'w', **kwargs) as dst:for i inrange(1, src.count + 1): source=rasterio.band(src, i), destination=rasterio.band(dst, i), src_transform=src.transform, resampling=Resampling.nearest
四、单城市演示计算（北京市，2020 年）
4.1 读取城市边界与人口栅格
city_all = gpd.read_file(PATH['city_shp'])city_all = city_all.to_crs(MYCRS)city_all['city_code'] = city_all['市代码']city_all['city_name'] = city_all['市']demo_cities = [110000, 120000, 130100]city_demo = city_all[city_all['city_code'].isin(demo_cities)]city_beijing = city_demo[city_demo['city_code'] == 110000]4.2 栅格转点函数
defraster_to_points(raster_path, city_geom):"""读取栅格数据，按城市边界裁剪，并转换为点数据 使用 all_touched=True 保留所有与边界多边形接触的像元， 与 R terra::mask(touches=TRUE) 的意图一致。 int32 栅格的 nodata 值无法被 np.isnan() 检测， 需先转为 float64 再替换 nodata 为 NaN。# all_touched=True：保留所有与边界多边形接触的像元# 与 R terra::mask(touches=TRUE) 的意图一致with rasterio.open(raster_path) as src: nodata_val = src.nodata # 通常是 -2147483648 (int32 最小值) out_image, out_transform = mask(src, [city_geom], crop=True, all_touched=True)# 将栅格数据转为 float64 以支持 NaN# 关键：int32 类型的栅格 nodata=-2147483648，np.isnan() 对 int 无效 data = out_image[0].astype(np.float64) data[data == nodata_val] = np.nan# 获取非空像元的坐标和值（只保留非 NaN 的格点） rows, cols = np.where(~np.isnan(data)) values = data[rows, cols] xs, ys = rasterio.transform.xy(out_transform, rows, cols) geometry=gpd.points_from_xy(xs, ys),# 只去掉 NaN，保留 pop=0 的格点（与 R 的 drop_na(pop) 一致）# pop=0 的格点会影响 LISA 分类中的中位数计算 gdf = gdf.dropna(subset=['pop']).reset_index(drop=True)pop_file = os.path.join(PATH['pop_tif'], "2020.tif")pts = raster_to_points(pop_file, city_beijing.geometry.iloc[0])4.3 构建空间权重矩阵
defbuild_spatial_weights(coords, threshold=1000): w = DistanceBand.from_array(coords, threshold=threshold, silence_warnings=True)coords = np.column_stack([pts.geometry.x, pts.geometry.y])w = build_spatial_weights(coords, ANALYSIS_PARAMS['dist_nb'])参数说明
transform = 'r'：行标准化权重。若一个格点有 4 个邻居，则每个邻居权重 = 1/4；有 2 个邻居则权重 = 1/2。
4.4 局部 Moran's I 检验
局部 Moran's I 的核心作用：找出人口高密度且周围也是高密度的区域（HH 聚集区），这正是城市中心的候选位置。
def_local_moran_conditional_p(y, w): 使用 R spdep::localmoran(conditional=TRUE) 的条件正态近似法计算 LISA p 值。 完全复现 R spdep::localmoran 源码中 conditional=TRUE 的计算逻辑： Ii = (zi / m2) * lag(zi) 其中 zi = y - mean(y), m2 = sum(z^2)/n E[Ii] = -(zi^2 * Wi) / ((n-1) * m2) 条件期望（与 unconditional 的 -Wi/(n-1) 不同） Var[Ii] = (zi/m2)^2 * (n/(n-2)) * (Wi2 - Wi^2/(n-1)) * (m2 - zi^2/(n-1)) z_score = (Ii - E[Ii]) / sqrt(Var[Ii]) p_i = 2 * (1 - Φ(|z_score|)) y = np.asarray(y, dtype=np.float64) m2 = np.sum(z ** 2) / n # mlvar=TRUE 时的方差 z_lag = lag_spatial(w, z)iflen(w.neighbors[i]) > 0: weights_i = np.array(w.weights[i]) Wi[i] = np.sum(weights_i) Wi2[i] = np.sum(weights_i ** 2)# 条件期望：E[Ii] = -(z_i^2 * Wi) / ((n-1) * m2) E_Ii = -(z ** 2 * Wi) / ((n - 1) * m2) var_Ii = (zi_over_m2 ** 2) * (n / (n - 2)) * (Wi2 - Wi ** 2 / (n - 1)) * (m2 - z ** 2 / (n - 1)) var_Ii = np.maximum(var_Ii, 0.0) z_score = np.where(var_Ii > 0, (Ii - E_Ii) / np.sqrt(var_Ii), 0.0) p_val = 2.0 * norm.sf(np.abs(z_score))deflocal_moran_lisa(pop_values, w, sig_level=0.05):"""计算局部 Moran's I 并进行 LISA 分类（条件正态近似法）""" p_values, Ii = _local_moran_conditional_p(pop_values, w) lag_val = lag_spatial(w, pop_values) med = np.median(pop_values) self_hl = np.where(pop_values > med, 'H', 'L') nbr_hl = np.where(lag_val > med, 'H', 'L') sig = p_values < sig_level types = np.where((sig) & (self_hl == 'H') & (nbr_hl == 'H'), 'HH', 'other')return p_values, lag_val, typesp_values, lag_val, types = local_moran_lisa(pts['pop'].values, w)pts['sig'] = p_values < ANALYSIS_PARAMS['sig_level']局部莫兰指数结果包含：
4.5 提取 HH 格点并聚类成中心
hh = pts[pts['type'] == 'HH'].copy()print(f"HH 格点数量: {len(hh)}")# 对 HH 格点再次做空间邻接，将连续区域标为同一 clusterdefcluster_centers(hh_gdf, dist_nb=1000):"""对 HH 格点进行空间聚类，识别城市中心""" coords = np.column_stack([hh_gdf.geometry.x, hh_gdf.geometry.y]) w_hh = DistanceBand.from_array(coords, threshold=dist_nb, silence_warnings=True) components = w_hh.component_labels hh_gdf['cluster'] = components centers = hh_gdf.groupby('cluster').agg({'geometry': lambda x: (np.mean(x.x), np.mean(x.y)) centers['n'] = hh_gdf.groupby('cluster').size().values centers['x'] = centers['geometry'].apply(lambda p: p[0]) centers['y'] = centers['geometry'].apply(lambda p: p[1]) centers = centers.drop('geometry', axis=1)centers = cluster_centers(hh, ANALYSIS_PARAMS['dist_nb'])centers = centers[(centers['n'] >= ANALYSIS_PARAMS['min_cells']) & (centers['pop'] >= ANALYSIS_PARAMS['min_pop'])]4.6 计算 5 大指标
五、批量串行计算（全国所有城市 × 所有年份）
5.1 封装核心计算函数
将以上步骤封装为函数，接收城市代码和人口 tif 文件路径，返回该城市该年份的 5 大指标：
defcalc_paper_indicators(city_code, pop_file, city_all, lw=None, params=None): city = city_all[city_all['city_code'] == city_code] city_geom = city.geometry.iloc[0] city_name = city['city_name'].iloc[0] pts = raster_to_points(pop_file, city_geom) coords = np.column_stack([pts.geometry.x, pts.geometry.y]) w = build_spatial_weights(coords, params['dist_nb'])# 局部 Moran's I 和 LISA 分类（条件正态近似法） p_values, lag_val, types = local_moran_lisa(pts['pop'].values, w, params['sig_level']) hh = pts[pts['type'] == 'HH'].copy()iflen(hh) < params['min_cells']: centers = cluster_centers(hh, params['dist_nb']) centers = centers[(centers['n'] >= params['min_cells']) & (centers['pop'] >= params['min_pop'])] pareto = calc_pareto_index(centers) poly = calc_polycentric_index(centers) cbd_idx = pts['pop'].idxmax() cbd_geom = pts.loc[cbd_idx, 'geometry'] sub_results = calc_decentralization(pts, cbd_geom) year_match = re.search(r'(\d{4})', os.path.basename(pop_file)) year = int(year_match.group(1)) if year_match elseNone'city_code': [city_code],'city_name': [city_name],'pareto': [round(pareto, 4)],'poly': [round(poly, 4)],'sub3': [round(sub_results['sub3'], 4)],'sub5': [round(sub_results['sub5'], 4)],'total_pop': [pts['pop'].sum()]print(f"计算失败: {city_code}, {pop_file}, 错误: {e}")result = calc_paper_indicators(110000, pop_file, city_demo)5.2 加载全部城市与年份文件
city_all = gpd.read_file(PATH['city_shp'])city_all = city_all.to_crs(MYCRS)city_all['city_code'] = city_all['市代码']city_all['city_name'] = city_all['市']pop_files = sorted(Path(PATH['pop_tif']).glob("*.tif"))print(f"找到 {len(pop_files)} 个年份文件")5.3 串行计算所有城市
六、改进算法：预计算空间权重矩阵
6.1 为什么可以改进？
在上面的批量计算中，对同一个城市，每个年份都重新计算了一次空间权重矩阵。但实际上，空间权重矩阵只取决于城市边界的形状——它与年份无关，只要城市行政边界不变（本项目使用 2021 年固定边界），同一城市所有年份的权重矩阵完全相同。
因此，可以先把所有城市的权重矩阵计算并保存好，然后在计算各年数据时直接读取，显著减少重复计算。
对于一个有 NNN 个城市、TTT 个年份的数据集：
当 T=25T = 25T=25（2000～2024 年）时，改进方法可将权重矩阵的计算量缩减为原来的 1/25。
6.2 预计算并保存所有城市的权重矩阵
defcompute_lw(city_code, city_all, sample_pop_file, output_dir="lwres"): city = city_all[city_all['city_code'] == city_code] city_geom = city.geometry.iloc[0] pts = raster_to_points(sample_pop_file, city_geom) coords = np.column_stack([pts.geometry.x, pts.geometry.y]) w = build_spatial_weights(coords, ANALYSIS_PARAMS['dist_nb']) os.makedirs(output_dir, exist_ok=True) output_path = os.path.join(output_dir, f"{city_code}.pkl")withopen(output_path, 'wb') as f:print(f"权重矩阵计算失败: {city_code}, 错误: {e}")# （joblib.Parallel 在 reticulate 环境下无法序列化 __main__ 中定义的函数，# 因此 Rmd 中使用串行循环；如需并行请使用独立的 batch_calculate.py 脚本）sample_pop_file = os.path.join(PATH['pop_tif'], "2020.tif")for code in city_demo['city_code'].unique(): compute_lw(code, city_demo, sample_pop_file)print("======== 所有城市权重矩阵计算完成，保存在 lwres/ ========")6.3 使用预计算权重矩阵的改进版计算函数
改进版函数从外部接收 lw 参数，不在函数内部重新计算权重矩阵：
defload_lw(city_code, lw_dir="lwres"): lw_path = os.path.join(lw_dir, f"{city_code}.pkl")if os.path.exists(lw_path):withopen(lw_path, 'rb') as f:defcalc_paper_indicators2(current_city, pop_file, lw, params=None): city_geom = current_city.geometry.iloc[0] city_code = current_city['city_code'].iloc[0] city_name = current_city['city_name'].iloc[0] pts = raster_to_points(pop_file, city_geom) p_values, lag_val, types = local_moran_lisa(pts['pop'].values, lw, params['sig_level']) hh = pts[pts['type'] == 'HH'].copy()iflen(hh) < params['min_cells']: centers = cluster_centers(hh, params['dist_nb']) centers = centers[(centers['n'] >= params['min_cells']) & (centers['pop'] >= params['min_pop'])] pareto = calc_pareto_index(centers) poly = calc_polycentric_index(centers) cbd_idx = pts['pop'].idxmax() cbd_geom = pts.loc[cbd_idx, 'geometry'] sub_results = calc_decentralization(pts, cbd_geom) year_match = re.search(r'(\d{4})', os.path.basename(pop_file)) year = int(year_match.group(1)) if year_match elseNone'city_code': [city_code],'city_name': [city_name],'pareto': [round(pareto, 4)],'poly': [round(poly, 4)],'sub3': [round(sub_results['sub3'], 4)],'sub5': [round(sub_results['sub5'], 4)],'total_pop': [pts['pop'].sum()]city_demo_single = city_demo[city_demo['city_code'] == 110000]lw_demo = load_lw(110000)result = calc_paper_indicators2(city_demo_single, pop_file, lw=lw_demo)6.4 使用改进算法批量计算全部数据
6.5 合并所有结果并保存
defmerge_results(input_dir="resb", output_file="results.dta"): all_files = list(Path(input_dir).glob("*.csv")) dfs = [pd.read_csv(f) for f in all_files] result = pd.concat(dfs, ignore_index=True)# 保存为 dta 格式（需要 pyreadstat） pyreadstat.write_dta(result, output_file)print(f"结果已保存为: {output_file}")# 如果没有 pyreadstat，保存为 CSV result.to_csv(output_file.replace('.dta', '.csv'), index=False)final_result = merge_results("resb", "2000～2024年各城市多中心指标（商玉萍版本）.dta")print(final_result.head())最终结果数据集包含如下变量：
| |
|---|
year | |
city_code | |
city_name | |
center | |
pareto | |
poly | |
sub3 | |
sub5 | |
total_pop | |
通常，pareto、poly 值越小，代表各中心的人口分布越均衡（均等）；sub3、sub5 越大，代表城市人口越去中心化。
七、稳健性检验：放宽中心人口门槛至 1 万人
商玉萍（2022）论文中提供了一项稳健性检验：将"总人口在 10 万人以上"的条件改为"总人口在 1 万人以上"，重新确定每个城市的中心数量（2center）作为稳健性指标。
这一操作的实现方式非常简单，只需将 ANALYSIS_PARAMS['min_pop'] 改为 10000，其余代码完全不变：
如何参加课程？
购买 RStata 名师讲堂会员即可参加该课程啦（之前的和未来的都可以参加）！
价格：2800/年 或者 4800/长期
购买会员可以从这里下单：https://rstata.duanshu.com/#/card/list/
名师讲堂会员权益：
- 参加平台上的其他 R 语言和 Stata 的课程；
- 以会员折扣价购买我们分享的数据资料（10 元/份）；
* 如果发票可添加小编微信 r_stata2 （RStata 李老师）开具。如需数据资料，购买后可添加小编微信免费领取数据折扣卡。
更多关于 RStata 会员的更多信息可添加微信号 r_stata2 咨询：
课程主页（点击文末的阅读原文即可跳转）：https://rstata.duanshu.com/#/brief/course/eb9f29edaadb45728360ba8f6aa0c399
