2010～2023 年各省市区县平均动态生境指数面板数据

2216 字 约 7 分钟 数据资料 面板数据生态环境

今天给大家分享一份 2010～2023 年各省市区县平均动态生境指数面板数据。该数据基于中国 500m 分辨率动态生境指数（DHI）栅格数据集处理得到，使用 R 语言的 terra 包按 2021 年行政区划（省 / 市 / 县）分区统计汇总为三级面板数据。

数据来源与处理方法

动态生境指数（Dynamic Habitat Indices, DHIs）从可用能量、环境稳定性、环境压力三个方面刻画生境生产力的年内物候特征，与生物多样性密切相关。本数据使用的原始栅格数据来自段昊玮等（2026）发表在《生态学报》上的《中国陆地生态系统动态生境指数格局及时空动态分析》一文所构建的数据集，其处理流程为：

  1. 基于 Google Earth Engine 平台获取 PML_V2_v0.1.8 蒸散发产品反演的 2010—2023 年全国总初级生产力（GPP）数据集，空间分辨率 500m，Albers 等积投影；
  2. 逐 8 天中位数合成月度 GPP，再计算逐年的三个动态生境指数（已排除沙漠、永久冰雪覆盖、水体等无生产力区域）；
  3. 该文的研究结论显示：空间上我国东南沿海向西北内陆累积生产力逐步递减、环境压力增加，东南沿海环境稳定性高于西北内陆；14 年间全国 58.49% 的地区 DHIcum 呈上升趋势。

在上述栅格数据的基础上，本数据所做的处理为分区域汇总：以 2021 年行政区划（省 / 市 / 县三级）为面要素，提取落入面要素范围内全部有效像元的算术平均值。需要注意：

  • 省 / 市 / 县三级是各自独立提取的（分别用对应层级的面要素直接统计），因此市级、省级结果是真实的像元加权均值，而不是由县级结果简单平均得到；
  • 统计口径为像元中心点落入面要素（exact = FALSE），区域面积过小（如部分市辖区）时可能没有任何像元中心点落入，结果记为缺失值 NA；
  • 原始栅格中被排除的沙漠、永久冰雪、水体等无生产力区域同样为缺失值。

数据概览

为了方便大家使用，我把栅格数据汇总成了省、市、县三级的 2010～2023 年平衡面板数据，包含 3 个指标：DHIcum（累计动态生境指数，年内 GPP 累积量）、DHImin（最小动态生境指数，年内 GPP 最小值）、DHIvar（动态生境指数变异度，年内 GPP 变异程度）。

三个文件的规模如下：

文件 观测值行数 区域单元数 时间范围 缺失观测数
2010～2023 年各省份平均动态生境指数.dta 476 34 个省份 2010～2023 年 0
2010～2023 年各城市平均动态生境指数.dta 5,194 371 个城市 2010～2023 年 14
2010～2023 年各区县平均动态生境指数.dta 40,278 2,877 个区县 2010～2023 年 238

变量说明：

  • 县级文件：省 省代码 市 市代码 县 县代码 年份 DHIcum DHImin DHIvar
  • 市级文件：省 省代码 市 市代码 年份 DHIcum DHImin DHIvar
  • 省级文件：省 省代码 年份 DHIcum DHImin DHIvar

缺失值含义：该年份该区域内没有落入任何有效像元（无生产力区域，或区域过小无像元中心点落入），例如北京市东城区这类面积较小的中心城区各年份均为缺失。

在用途上，DHIs 在生态环境质量监测、生物多样性管理与评估中具有不同的侧重和互补的重要意义；这类县域面板数据也可直接用于区域层面的实证分析，例如参考文献中以中国光伏扶贫项目为准自然实验、基于 2010—2024 年全国县域面板数据采用多时点双重差分模型评估光伏项目建设对县域生物多样性影响的研究。

数据预览如下（依次为区县、城市、省份三级）：

图表展示

下图展示了 2023 年各省份累计动态生境指数（DHIcum）的分布：

下图展示了 2023 年各区县动态生境指数变异度（DHIvar）的分布：

下图展示了 2010—2023 年三个 DHI 指标的时间趋势（灰色细线为各省份，红色粗线为全国均值，蓝色虚线为全国均值的 OLS 线性趋势，左上角标注了年均变化量）：

下图展示了 2023 年县级 DHIcum 与 DHIvar 的关系（散点按 DHImin 着色，并给出 OLS 拟合线、R² 与 Pearson 相关系数）：

处理代码

汇总过程使用 R 语言的 terra、sf、tidyverse、haven 包完成：先把同一年份的 3 个指标合并为一个三波段 SpatRaster 一次读入，再对省 / 市 / 县三组面要素分别做分区统计，最后加上变量标签写出 Stata 15 格式的 dta 文件。关键代码如下：

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

# 逐年分区统计：三个指标合成一个三波段栅格，一次读入同时提取
for (y in years) {
  r <- rast(sapply(dhi_vars, dhi_path, y = y))   # DHIcum / DHImin / DHIvar
  names(r) <- dhi_vars

  for (lv in names(levels_def)) {
    ld <- levels_def[[lv]]
    # 以像元中心点落入面要素统计区域内有效像元的算术均值
    e  <- terra::extract(r, ld$v, fun = mean, na.rm = TRUE, ID = FALSE,
                         exact = FALSE)
    if (nrow(e) != nrow(ld$key)) stop("提取结果与要素数量不一致：", lv, " ", y)

    res_list[[length(res_list) + 1]] <- bind_cols(ld$key, as_tibble(e)) |>
      mutate(年份 = y, 层级 = lv)
  }
  rm(r); invisible(gc(verbose = FALSE))
}

# NaN 统一转为 NA，然后拆分为省 / 市 / 县三个数据集
res_all <- bind_rows(res_list) |>
  mutate(across(all_of(dhi_vars), ~ ifelse(is.nan(.x), NA_real_, .x)),
         年份 = as.integer(年份))

# 加变量标签并写出 dta（Stata 15 格式）
write_dta(add_var_labels(df_cnty), f_cnty, version = 15,
          label = "数据处理：微信公众号 RStata")

图表绘制使用 ggplot2 + sf 完成，地图底图来自 2021 年行政区划 mini 版，配色使用 ColorBrewer 的 YlGnBu / YlOrRd 色带，并添加了比例尺与指北针：

# 按分位数分组 + 生成配色（含“无数据”灰色组）
make_breaks <- function(x, n = 6, digits = 3) {
  qs   <- unique(quantile(x, probs = seq(0, 1, length.out = n + 1), na.rm = TRUE))
  n2   <- length(qs) - 1
  labs <- sprintf("%s ~ %s",
                  sprintf("%.3f", qs[1:n2]), sprintf("%.3f", qs[2:(n2 + 1)]))
  g <- as.character(cut(x, breaks = qs, include.lowest = TRUE, labels = labs))
  g[is.na(g) & is.na(x)] <- "无数据"     # NA 归为“无数据”组，配色为灰色
  factor(g, levels = c(labs, if (any(is.na(x))) "无数据"))
}

ggplot(provmap_data) +
  geom_sf(data = plotbbox, fill = "#BBD1EB", color = NA) +
  geom_sf(aes(fill = DHIcum_g), color = "grey30", linewidth = 0.15) +
  geom_sf(data = china_neighboring, fill = "#EDEDED", color = "grey70", linewidth = 0.1) +
  scale_fill_manual(name = "2023 年 DHIcum", values = pal1, na.value = "grey85",
                    drop = FALSE) +
  annotation_scale(location = "bl", width_hint = 0.25, text_family = cnfont) +
  annotation_north_arrow(location = "tr", which_north = "false",
                         style = north_arrow_fancy_orienteering(text_family = cnfont)) +
  labs(title = "2023 年中国各省份累计动态生境指数（DHIcum）",
       subtitle = "数据处理&绘制：微信公众号 RStata",
       caption = "数据来源：RStata 数据中心: 2010～2023 年各省市区县平均动态生境指数面板数据. 2026. https://tidyfriday.cn/rsdb2/") -> p1

ggsave("picdir/图1_2023年各省份累计动态生境指数DHIcum分布.png",
       plot = p1, width = 10, height = 8, device = png, bg = "white")

附件中也提供了该数据的处理代码供参考：

数据引用格式

由于该数据包含较多 RStata 处理的内容，在研究中使用该数据请使用清晰的方式注明数据来源于 RStata 或者 RStata 数据中心，并使用如下格式引用：

RStata 数据中心: 2010～2023 年各省市区县平均动态生境指数面板数据. 2026. https://tidyfriday.cn/rsdb2/

英文文献可以使用下面的格式引用：

RStata Data Center: Panel Data on Average Dynamic Habitat Indices (DHI) across Provinces, Cities, and Districts/Counties in China, 2010–2023. 2026. https://tidyfriday.cn/rsdb2/

原始栅格数据集的参考文献：

段昊玮, 杨莹莹, 万华伟, 等. 中国陆地生态系统动态生境指数格局及时空动态分析[J]. 生态学报, 2026, 46(2): 691-704. DOI: 10.20103/j.stxb.202505191229

关联课程/数据推荐

1985～2025 年各城市核心绿色生境景观破碎化指标面板数据：https://rstata.duanshu.com/#/course/bdfd89c258af4b2bad8d9e73b75d2176

2001～2023 年中国各省市区县生态环境质量面板数据：https://rstata.duanshu.com/#/course/865c987f170c4bc48debe58c610fdd7a

2000 年 2 月～2024 年 12 月各省市区县平均植被覆盖度：https://rstata.duanshu.com/#/course/8da9d36689d6465896ce298a77c34156

RStata 定制｜栅格数据转面板数据或时间序列数据：https://rstata.duanshu.com/#/course/839afebb01d54457b42b94e87e7dc085

如何将栅格数据处理成面板数据或时间序列数据？以 PM2.5 浓度数据处理为例（含并行计算内容）：https://rstata.duanshu.com/#/course/109dfcdfd1174d1eb2d00c621df4f047

如果有相关需要可以联系李老师付费定制。

点击这里跳转到 RStata 短书平台获取附件：2010～2023 年各省市区县平均动态生境指数面板数据

相关文章

评论