2000～2023 年上市公司周边一定范围的 NDVI 均值和最大值面板数据

今天给大家分享一份 2000～2023 年上市公司周边一定范围的 NDVI 均值和最大值面板数据。原始 NDVI 栅格数据来源于 MODIS 卫星遥感数据：

2000～2023 年各省市区县 NDVI 面板数据（使用月度最大值和均值合成）：https://rstata.duanshu.com/#/course/0342f29770b74d639c2a73e241c3c830

NDVI（归一化植被指数，Normalized Difference Vegetation Index）是反映植被生长状况的重要遥感指标，取值范围为 -1 到 1，值越高表示植被覆盖度越高。原始栅格数据分辨率为 250m，由于直接基于 250m 栅格计算非常耗时，因此先将栅格数据聚合到 1km 分辨率再进行计算。

上市公司地址数据使用了之前分享的：

2000～2024 年上市公司注册地址与办公地址（含经纬度及其所处的省市区县）：https://rstata.duanshu.com/#/course/746d4f595eba41a2a116b5a6edc3ac0a

数据概览

为便于大家使用，我计算了 2000～2023 年上市公司办公地址和注册地址周边 1km、2km、3km、4km、5km、10km、15km、20km、25km、30km 范围内的平均 NDVI 和最大 NDVI，最终得到了附件中的面板数据。

数据包含如下变量：

股票代码
股票简称
年份
地址类型（办公地址 / 注册地址）
距离范围（1、2、3、4、5、10、15、20、25、30）
平均NDVI（缓冲区范围内 NDVI 均值，原始值 ×10000）
最大NDVI（缓冲区范围内 NDVI 最大值，原始值 ×10000）

经过计算之后一共得到了 126 万条结果，覆盖 2000～2023 年共 24 个年份。

数据预览如下：

附件中还包含了上市公司地址数据：

图表展示

下图展示了 2000～2023 年上市公司周边 5km 范围内平均 NDVI 的变化趋势（含 95% 置信区间）：

下图展示了 2020 年上市公司办公地址周边 5km 范围内平均 NDVI 与最大 NDVI 的关系：

下图展示了 2020 年各省份上市公司周边平均 NDVI 的空间分布：

处理代码

该数据使用 R 语言的 terra、sf 包进行处理，通过 Stata 绘制展示图。处理流程如下：

第一步：将原始 250m 栅格数据重采样到 1km 分辨率

原始栅格数据分辨率约为 250m（0.002499889°），聚合 factor = 4 可得到约 1km 分辨率。其中月度均值合成栅格用 mean 聚合，月度最大值合成栅格用 max 聚合：

# 原始分辨率 0.002499889°（≈278 m），聚合 factor = 4 → ~1 km
# rawdata3：月度均值合成 → 用 mean 聚合
rst <- terra::rast(f3)
terra::aggregate(rst,
                 fact  = 4,
                 fun   = "mean",
                 na.rm = TRUE) %>%
  terra::writeRaster(o3, overwrite = TRUE)

# rawdata2：月度最大值合成 → 用 max 聚合
rst <- terra::rast(f2)
terra::aggregate(rst,
                 fact  = 4,
                 fun   = "max",
                 na.rm = TRUE) %>%
  terra::writeRaster(o2, overwrite = TRUE)

第二步：基于 1km 栅格计算缓冲区内 NDVI

使用并行计算（5 核），对每个上市公司地址的缓冲区分别提取平均 NDVI 和最大 NDVI：

# 并行计算：对每个缓冲区文件分别提取平均 NDVI 和最大 NDVI
parLapply(cl, 1:nrow(yddf), function(x) {
  readr::read_rds(yddf$value[x]) -> dfbuffer1
  dfbuffer1 %>%
    sf::st_drop_geometry() %>%
    dplyr::mutate(ID = dplyr::row_number()) -> temp1

  yr <- yddf$year[x]

  # 平均 NDVI：使用 1 km 均值合成栅格（rawdata3_1km）
  terra::rast(file.path("rawdata3_1km", paste0(yr, ".tif"))) -> rst_mean
  terra::extract(rst_mean, terra::vect(dfbuffer1),
                 fun = "mean", na.rm = TRUE, exact = F) %>%
    dplyr::as_tibble() %>%
    dplyr::select(-ID) %>%
    dplyr::rename(平均NDVI = 1) -> mean_val

  # 最大 NDVI：使用 1 km 最大值合成栅格（rawdata2_1km）
  terra::rast(file.path("rawdata2_1km", paste0(yr, ".tif"))) -> rst_max
  terra::extract(rst_max, terra::vect(dfbuffer1),
                 fun = "max", na.rm = TRUE, exact = F) %>%
    dplyr::as_tibble() %>%
    dplyr::select(-ID) %>%
    dplyr::rename(最大NDVI = 1) -> max_val

  # 合并并保存
  bind_cols(temp1, mean_val, max_val) %>%
    dplyr::mutate(file = yddf$value[x]) %>%
    readr::write_rds(paste0("res/", x, ".rds"))
})

第三步：合并所有分块结果并整理列结构

# 合并所有分块结果
parLapply(cl, ls_res, read_rds) %>%
  bind_rows() -> df

# 整理列结构
df %>%
  select(股票代码, file, 平均NDVI, 最大NDVI) %>%
  mutate(file = basename(file)) %>%
  tidyr::extract(
    file,
    into    = c("year", "type", "dist"),
    regex   = "(.*)_(.*)_(.*)\\.rds",
    convert = TRUE,
    remove  = TRUE
  ) -> df2

# 补充股票简称
df2 %>%
  rename(年份 = year) %>%
  left_join(dfo) %>%
  rename(股票简称1 = 股票简称) %>%
  left_join(dfo2) %>%
  mutate(股票简称1 = if_else(is.na(股票简称1), 股票简称, 股票简称1)) %>%
  select(-股票简称) %>%
  rename(股票简称 = 股票简称1, 地址类型 = type, 距离范围_km = dist) -> df3

# 保存为 dta
df3 %>%
  haven::write_dta(
    "2000～2023年上市公司周边一定范围的 NDVI 均值和最大值面板数据.dta",
    label = "数据处理：微信公众号 RStata"
  )

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

数据引用格式

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

RStata 数据中心: 2000～2023年上市公司周边一定范围的 NDVI 均值和最大值面板数据. 2026. https://tidyfriday.cn/rsdb2/

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

RStata Data Center: Panel Data on Mean and Maximum NDVI within Certain Ranges around Listed Companies, 2000–2023. 2026. https://tidyfriday.cn/rsdb2/

关联课程/数据推荐

2000～2023 年各省市区县 NDVI 面板数据（使用月度最大值和均值合成）：https://rstata.duanshu.com/#/course/0342f29770b74d639c2a73e241c3c830

2000～2024 年上市公司注册地址与办公地址（含经纬度及其所处的省市区县）：https://rstata.duanshu.com/#/course/746d4f595eba41a2a116b5a6edc3ac0a

使用 R 语言计算工企周边一定范围内的平均 PM2.5 浓度：https://rstata.duanshu.com/#/course/016330182fb64ecd9fab8ada0d816130

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

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

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

点击这里跳转到 RStata 短书平台获取附件：2000～2023 年上市公司周边一定范围的 NDVI 均值和最大值面板数据

数据概览

图表展示

处理代码

数据引用格式

关联课程/数据推荐

相关文章

评论