名师讲堂｜使用 Python 测算城乡夜间灯光亮度及城乡差异

3025 字 约 10 分钟 名师讲堂 夜间灯光Python

今天给大家分享使用 Python 测算中国各地区城乡夜间灯光亮度差异的方法。该方法基于城市建成区与镇建成区的矢量边界数据，对夜间灯光遥感栅格数据进行分区统计，计算城市与农村的平均灯光亮度及其比值。

原理与方法

城乡夜间灯光亮度差异是衡量城镇化水平和城乡差距的重要指标。本方法的核心思想是：

  1. 将城市建成区（Cities）和镇建成区（Towns）的矢量边界合并作为”城市”区域
  2. 其余区域定义为”农村”区域
  3. 对夜间灯光栅格数据进行分区统计，计算城市和农村的平均亮度
  4. 用城乡亮度比值（城市/农村）衡量城乡差距

其中城市建成区和镇建成区的数据来源于这个论文：

Bai, M., Zhang, X., Ai, W., Jing, X., & Liu, L. (2026). A high-resolution global annual city and town boundaries dataset (2000–2022) derived from GLC_FCS30D product. Scientific Data, 13:50. https://doi.org/10.1038/s41597-025-02568-9

夜间灯光数据平台上有好几个版本，这里选择的是施开放版本：

1992～2024 年中国各省市区县乡镇类 DMSP-OLS 夜间灯光面板数据（施开放版本）：https://rstata.duanshu.com/#/brief/course/2c3880893f5947a5bf764befd3994c74

另外也可以选择其余的三个版本：

2000～2024 年各省市区县乡镇类 NPP-VIIRS 夜间灯光面板数据 & 栅格数据（余柏蒗版本1）：https://rstata.duanshu.com/#/brief/course/f8e6c2d0d4e549a58749d7d86a60e238

1992～2024 年各省市区县乡镇类 NPP-VIIRS 夜间灯光面板数据 & 栅格数据（余柏蒗版本2）：https://rstata.duanshu.com/#/brief/course/5ee752e4c1d84ab48b2235366dfb83ec

1986～2024 年中国各省市区县、乡镇夜间灯光面板数据（田一禾版本）：https://rstata.duanshu.com/#/brief/course/b7392dc0a7f649c9bf129666445c866e

指标定义

指标 定义 说明
城市平均亮度 建成区栅格的灯光均值 值越大表示城镇化水平越高
农村平均亮度 非建成区栅格的灯光均值 值越大表示农村发展水平越高
城乡比 城市/农村亮度比值 值越大表示城乡差距越大

技术流程

栅格数据 (TIF)
       ↓
  投影转换 (EPSG:4326)
       ↓
  建成区掩膜 (Cities + Towns)
       ↓
  城市灯光 × 掩膜
  农村灯光 × (1 - 掩膜)
       ↓
  分区统计 (zonal mean)
       ↓
  省 / 市 / 区县 面板数据

使用 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_light"
.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} 代码块时，reticulate 已经通过 RETICULATE_PYTHON 环境变量知道要使用 .venv_light，不会再去碰 Anaconda。

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

# 检查关键包是否已安装，缺失的才安装
py_pkgs <- c(
  "numpy", "pandas", "geopandas", "rasterio", "rasterstats"
)
installed <- py_list_packages(".venv_light")$package
need_install <- setdiff(py_pkgs, installed)

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

验证激活状态

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

查看已安装的包

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

虚拟环境管理常用命令

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

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

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

数据准备

本方法需要以下数据：

  1. 夜间灯光栅格数据：年度 TIF 格式的夜间灯光亮度数据（以施开放版本为例）
  2. 城市建成区边界：Cities_2022.shp（城市建成区矢量数据）
  3. 镇建成区边界：Towns_2022.shp（镇建成区矢量数据）
  4. 行政区划边界：省/市/区县 shp 文件（用于分区统计）

数据来源说明：

  • 夜间灯光栅格数据：可使用施开放版本、田一禾版本、余柏蒗版本等
  • 建成区边界数据：来自 Cities_2000_2022 和 Towns_2000_2022 文件夹
  • 行政区划数据：2021 年行政区划

首先设置参数：

import os
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings('ignore')

# 设置参数
VERSION_NAME = "施开放版本"
VERSION_FOLDER = "夜间灯光亮度栅格数据（施开放版本）"
YEAR = 2022

print("=" * 60)
print("城乡夜间灯光亮度差异测算 - Python 版本")
print(f"数据版本: {VERSION_NAME}")
print(f"年份: {YEAR}")
print("=" * 60)

读取灯光栅格数据

加载夜间灯光栅格数据并投影到 WGS84：

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

# 读取灯光栅格数据
light_tif_path = os.path.join(VERSION_FOLDER, f"{YEAR}.tif")
light_raster = rasterio.open(light_tif_path)

print(f"栅格文件: {light_tif_path}")
print(f"原始CRS: {light_raster.crs}")
print(f"原始维度: {light_raster.height} x {light_raster.width}")
print(f"数据范围: {light_raster.bounds}")

# 读取数据
light_data = light_raster.read(1)
print(f"数值范围: {np.nanmin(light_data)} - {np.nanmax(light_data)}")

投影到 WGS84（用于与 shp 数据匹配）：

import rasterio.transform

# 投影到 WGS84 (EPSG:4326)
light_crs_epsg = light_raster.crs.to_epsg()
print(f"原始 EPSG: {light_crs_epsg}")

if light_crs_epsg != 4326:
    print("正在投影到 WGS84 (EPSG:4326)...")

    dst_crs = 'EPSG:4326'
    transform, width, height = calculate_default_transform(
        light_raster.crs, dst_crs, light_raster.width, light_raster.height, *light_raster.bounds
    )

    light_proj_data = np.empty((height, width), dtype=light_data.dtype)

    reproject(
        source=light_data,
        destination=light_proj_data,
        src_transform=light_raster.transform,
        src_crs=light_raster.crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.bilinear
    )

    print(f"投影后维度: {height} x {width}")
else:
    light_proj_data = light_data
    transform = light_raster.transform
    print("数据已经是 WGS84 坐标系，无需投影")

构建城乡掩膜

城乡掩膜的构建是整个方法的核心步骤：

import geopandas as gpd
from rasterio.features import rasterize
from rasterio.transform import from_bounds

# 加载建成区边界数据
此处代码需下载讲义材料查看～

将矢量边界栅格化为掩膜：

from rasterio.coords import BoundingBox

# 获取栅格参数
# 此处代码需下载讲义材料查看～

print(f"掩膜城市像元数: {np.sum(mask_raster_data == 1)}")
print(f"掩膜农村像元数: {np.sum(mask_raster_data == 0)}")

关键说明：不再按面积过滤建成区，而是直接合并所有城市和镇矢量区域。这是因为 Towns 数据中存在面积字段为 0 的记录（可能是原始数据采集问题），去掉面积过滤可以保留更多有效区域。

计算城乡灯光

将灯光栅格与掩膜相乘，分别得到城市和农村的灯光强度：

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

分区统计

使用 rasterstats.zonal_stats 函数对省、市、县分别进行分区统计，计算各区划内的灯光均值。

读取行政区划数据

# 尝试读取行政区划数据
if (dir.exists("2021行政区划")) {
  province_shp <- sf::st_read("2021行政区划/省.shp", quiet = TRUE)
  city_shp <- sf::st_read("2021行政区划/市.shp", quiet = TRUE)
  county_shp <- sf::st_read("2021行政区划/县.shp", quiet = TRUE)

  cat("行政区划数据读取完成\n")
  cat("- 省:", nrow(province_shp), "个\n")
  cat("- 市:", nrow(city_shp), "个\n")
  cat("- 县:", nrow(county_shp), "个\n")
} else {
  cat("未找到行政区划数据，将使用建成区数据中的属性信息\n")
}

省级统计

from rasterstats import zonal_stats

# 此代码块仅作示例，实际运行需要行政区划数据
# 此处代码需下载讲义材料查看～

# 提取结果
prov_data_list = []
for feat_u, feat_r in zip(prov_urban_stats, prov_rural_stats):
    props_u = feat_u['properties']
    props_r = feat_r['properties']

    urban_mean = props_u.get('mean', np.nan)
    rural_mean = props_r.get('mean', np.nan)

    # 计算城乡比
    if rural_mean and rural_mean > 0 and urban_mean and not np.isnan(urban_mean):
        ratio = urban_mean / rural_mean
    else:
        ratio = np.nan

    prov_data_list.append({
        '年份': YEAR,
        '省代码': props_u.get('省代码', ''),
        '省': props_u.get('省', ''),
        '城市平均亮度': urban_mean if urban_mean is not None else np.nan,
        '农村平均亮度': rural_mean if rural_mean is not None else np.nan,
        '城乡比': ratio
    })

prov_data = pd.DataFrame(prov_data_list)
print(prov_data.head(10).to_string(index=False))

完整计算代码

以下是完整的 Python 计算代码，可以直接运行：

#!/usr/bin/env python3
# 完整代码请参考同目录下的 main.py 文件
# 运行命令: python main.py

运行 Python 脚本

在 R 中直接调用 Python 脚本进行计算：

# 运行 Python 脚本
py_run_file("main.py")

读取结果

# 读取计算结果
result_folder <- "施开放版本"
version_name <- "施开放版本"
year <- 2022

if (dir.exists(result_folder)) {
  # 读取省级结果
  prov_file <- file.path(result_folder, sprintf("%s省级结果_%d.csv", version_name, year))
  if (file.exists(prov_file)) {
    prov_data <- read.csv(prov_file, encoding = "UTF-8")
    cat("省级结果:\n")
    print(head(prov_data, 10))
  }

  # 读取市级结果
  city_file <- file.path(result_folder, sprintf("%s市级结果_%d.csv", version_name, year))
  if (file.exists(city_file)) {
    city_data <- read.csv(city_file, encoding = "UTF-8")
    cat("\n市级结果条数:", nrow(city_data), "\n")
  }

  # 读取区县级结果
  county_file <- file.path(result_folder, sprintf("%s区县级结果_%d.csv", version_name, year))
  if (file.exists(county_file)) {
    county_data <- read.csv(county_file, encoding = "UTF-8")
    cat("区县级结果条数:", nrow(county_data), "\n")
  }
}

结果描述性统计

查看各指标的描述性统计：

if (exists("prov_data")) {
  cat("=== 省级结果描述性统计 ===\n")
  print(summary(prov_data[, c("城市平均亮度", "农村平均亮度", "城乡比")]))
}

if (exists("city_data")) {
  cat("\n=== 市级结果描述性统计 ===\n")
  print(summary(city_data[, c("城市平均亮度", "农村平均亮度", "城乡比")]))
}

if (exists("county_data")) {
  cat("\n=== 区县级结果描述性统计 ===\n")
  print(summary(county_data[, c("城市平均亮度", "农村平均亮度", "城乡比")]))
}

结果变量说明

最终输出的面板数据包含以下变量：

变量名 说明 省级 市级 区县
年份 观测年份 ✓ ✓ ✓
省代码 省级区划代码 ✓ ✓ ✓
省 省级区划名称 ✓ ✓ ✓
市代码 市级区划代码 ✓ ✓
市 市级区划名称 ✓ ✓
县代码 县级区划代码 ✓
县 县级区划名称 ✓
城市平均亮度 建成区平均灯光强度 ✓ ✓ ✓
农村平均亮度 非建成区平均灯光强度 ✓ ✓ ✓
城乡比 城市/农村亮度比值 ✓ ✓ ✓

R 与 Python 函数对照

以下是本方法中主要函数的 R 与 Python 对照表：

R 函数 Python 函数 说明
terra::rast() rasterio.open() 读取栅格数据
terra::project() rasterio.warp.reproject() 投影转换
sf::st_read() geopandas.read_file() 读取矢量数据
sf::st_transform() gdf.to_crs() 坐标转换
sf::st_make_valid() gdf.make_valid() 修复几何
terra::rasterize() rasterio.features.rasterize() 栅格化
terra::zonal() rasterstats.zonal_stats() 分区统计
terra::values() numpy.array 提取栅格值
dplyr::left_join() pandas.merge(how='left') 左连接
haven::write_dta() pandas.to_csv() 保存数据

关键差异说明

1. 栅格数据处理

  • R (terra): 使用 SpatRaster 对象，支持延迟计算
  • Python (rasterio): 使用 NumPy 数组，直接操作内存数据

2. 分区统计

  • R (terra): zonal() 函数直接返回数据框
  • Python (rasterstats): zonal_stats() 返回字典列表，需要额外处理

3. 坐标参考系统

  • R: 使用 EPSG 代码字符串，如 “EPSG:4326”
  • Python: rasterio 使用 CRS 对象，geopandas 支持 EPSG 代码

4. 缺失值处理

  • R: 使用 NA
  • Python: 使用 numpy.nan

附录：完整 Python 代码

完整的 Python 计算代码已保存为 main.py，可以直接运行：

python main.py

或在 R 中调用：

reticulate::py_run_file("main.py")

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Python 测算城乡夜间灯光亮度及城乡差异

相关文章

评论