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

原理与方法

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

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

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

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)
       ↓
  省 / 市 / 区县 面板数据

数据准备

本方法需要以下数据：

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

数据来源说明：

夜间灯光栅格数据：可使用施开放版本、田一禾版本、余柏蒗版本等

建成区边界数据：来自 Cities_2000_2022 和 Towns_2000_2022 文件夹

行政区划数据：2021 年行政区划

首先加载必要的 R 包并设置参数：

library(tidyverse)
library(sf)
library(terra)
library(haven)

setwd("/Users/ac/Desktop/城乡夜间灯光强度对比")

# 以施开放版本 2022年为例
version_name <- "施开放版本"
version_folder <- "夜间灯光亮度栅格数据（施开放版本）"
year <- 2022

读取行政区划数据

首先读取行政区划矢量边界，用于后续分区统计：

province_shp <- st_read("2021行政区划/省.shp", quiet = TRUE)
city_shp <- st_read("2021行政区划/市.shp", quiet = TRUE)
county_shp <- st_read("2021行政区划/县.shp", quiet = TRUE)

cat("行政区划数据读取完成\n")
cat("- 省:", nrow(province_shp), "个\n")
cat("- 市:", nrow(city_shp), "个\n")
cat("- 县:", nrow(county_shp), "个\n")

# 提取行政区划参照表（用于后续补充变量）
# 市-省对应关系
city_province_df <- city_shp %>%
  st_drop_geometry() %>%
  select(市代码, 市, 省代码, 省) %>%
  distinct(市代码, .keep_all = TRUE)

# 县-市-省对应关系
county_city_province_df <- county_shp %>%
  st_drop_geometry() %>%
  select(县代码, 县, 市代码, 市, 省代码, 省) %>%
  distinct(县代码, .keep_all = TRUE)

cat("市级参照表:", nrow(city_province_df), "条\n")
cat("区县级参照表:", nrow(county_city_province_df), "条\n")

读取灯光栅格数据

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

# 加载灯光栅格数据
light_raster <- rast(file.path(version_folder, sprintf("%s.tif", year)))
cat("原始灯光栅格维度:", dim(light_raster), "\n")
cat("原始灯光范围:", paste(as.character(minmax(light_raster)), collapse = " - "), "\n")

# 投影到 WGS84（用于与 shp 数据匹配）
light_proj <- project(light_raster, "EPSG:4326", method = "bilinear")
cat("投影后栅格维度:", dim(light_proj), "\n")
cat("投影后范围:")
ext(light_proj)

构建城乡掩膜

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

# 加载建成区边界数据
cities <- st_read(sprintf("Cities_2000_2022/Cities_%s.shp", year), quiet = TRUE)
towns <- st_read(sprintf("Towns_2000_2022/Towns_%s.shp", year), quiet = TRUE)

# 合并所有城市和镇作为城市区域（不进行面积过滤）
urban <- rbind(
  cities %>% select(geometry),
  towns %>% select(geometry)
)

cat("城市建成区:", nrow(cities), "个\n")
cat("镇建成区:", nrow(towns), "个\n")
cat("合并后城市区域总数:", nrow(urban), "个\n")

将矢量边界转换为 terra 格式并栅格化为掩膜：

# 转为 terra 格式（修复无效几何）
# 此处代码需下载讲义材料查看～

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

计算城乡灯光

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

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

分区统计

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

辅助函数

首先定义一个辅助函数，用于提取 zonal 统计结果：

zonal_to_df <- function(z_result, code_col) {
  df <- as.data.frame(z_result)
  names(df)[1] <- code_col
  names(df)[2] <- "value"
  df
}

省级统计

以省为单位进行灯光均值统计：

# 省级分区栅格
province_proj <- st_transform(province_shp, crs = "EPSG:4326")
province_raster <- rasterize(vect(province_proj), light_proj, field = "省代码")

# 计算各省的城市/农村灯光均值
prov_urban <- zonal_to_df(zonal(urban_lights, province_raster, fun = "mean", na.rm = TRUE), "省代码")
prov_rural <- zonal_to_df(zonal(rural_lights, province_raster, fun = "mean", na.rm = TRUE), "省代码")

# 合并结果
prov_merge <- province_proj %>% st_drop_geometry() %>% select(省代码, 省)
prov_data <- prov_merge %>%
  left_join(rename(prov_urban, urban_mean = value), by = "省代码") %>%
  left_join(rename(prov_rural, rural_mean = value), by = "省代码")

# 添加年份和计算城乡比
prov_data$year <- as.integer(year)
prov_data$ratio <- prov_data$urban_mean / prov_data$rural_mean
# 处理异常值（农村亮度为0或负数时城乡比为NA）
prov_data$ratio[is.infinite(prov_data$ratio) | is.nan(prov_data$ratio) | prov_data$rural_mean <= 0] <- NA

# 整理列名
colnames(prov_data) <- c("省代码", "省", "城市平均亮度", "农村平均亮度", "年份", "城乡比")
prov_data <- prov_data %>% select(年份, 省代码, 省, 城市平均亮度, 农村平均亮度, 城乡比)

cat("\n省级结果示例（", year, "年）：\n")
prov_data %>% head(10)

市级统计

市级统计方法与省级相同：

# 市级分区栅格
# 此处代码需下载讲义材料查看～

区县级统计

区县级统计同样适用：

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

cat("区县结果条数:", nrow(county_data), "\n")
county_data %>% head(10)

保存结果

将计算结果保存为 Stata 格式：

# 创建结果文件夹
result_folder <- version_name
dir.create(result_folder, showWarnings = FALSE)

# 省级结果
write_dta(prov_data, file.path(result_folder, sprintf("%s省级结果_%d.dta", version_name, year)),
          version = 14, label = "数据处理：微信公众号 RStata")

# 市级结果
write_dta(city_data, file.path(result_folder, sprintf("%s市级结果_%d.dta", version_name, year)),
          version = 14, label = "数据处理：微信公众号 RStata")

# 区县级结果
write_dta(county_data, file.path(result_folder, sprintf("%s区县级结果_%d.dta", version_name, year)),
          version = 14, label = "数据处理：微信公众号 RStata")

cat("结果已保存到", result_folder, "文件夹\n")

结果描述性统计

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

cat("=== 省级结果描述性统计 ===\n")
prov_data %>% select(城市平均亮度, 农村平均亮度, 城乡比) %>% summary()

cat("\n=== 市级结果描述性统计 ===\n")
city_data %>% select(城市平均亮度, 农村平均亮度, 城乡比) %>% summary()

cat("\n=== 区县级结果描述性统计 ===\n")
county_data %>% select(城市平均亮度, 农村平均亮度, 城乡比) %>% summary()

结果变量说明

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

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

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