2001～2025 年上市公司周边的 NDVI 面板数据

之前给大家分享过 2001～2025 年各省市区县月度与年度NDVI面板数据，不少小伙伴希望把 NDVI 进一步精确到企业周边。今天给大家分享一份 2001～2025 年上市公司周边的 NDVI 面板数据，由 RStata 数据中心处理得到（基于上市公司地址经纬度与年度 NDVI 栅格提取，自制数据）。

2001～2025 年各省市区县月度与年度NDVI面板数据：https://rstata.duanshu.com/#/brief/course/a431c448b7294c8598604a006d37f3e2

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

数据概览

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

该数据共包含 1,465,030 条观测，覆盖 5,844 家上市公司、2001～2025 年、两种地址类型（办公地址 / 注册地址）以及 10 个距离圈层。数据概览：

包含如下变量：

股票代码、股票简称、年份、地址类型、距离_km、平均NDVI、最大NDVI

变量含义：

股票代码 / 股票简称：上市公司的标识信息；
年份：2001～2025 年；
地址类型：办公地址或注册地址（上市公司办公地与注册地的经纬度不同，因此分别构建缓冲区）；
距离_km：以公司地址为中心、向外缓冲的距离圈层半径，取值为 1、2、3、4、5、10、15、20、25、30（单位：km）；
平均NDVI：该距离圈层缓冲区范围内，由年度均值合成 NDVI 栅格提取得到的平均植被指数；
最大NDVI：该距离圈层缓冲区范围内，由年度最大值合成 NDVI 栅格提取得到的最大植被指数。

该数据可用于刻画上市公司办公 / 注册地址周边不同半径范围内的植被覆盖状况，为研究企业所在地的生态环境质量、城市绿化空间分布，及其与公司行为（如绿色创新、环境绩效、区位选择等）之间的关系，提供微观地理层面的度量。

图表展示

下图展示了 2001～2025 年上市公司周边平均 NDVI 与最大 NDVI 的逐年变化趋势：

下图展示了平均 NDVI 与最大 NDVI 的关系，并叠加了 OLS 拟合线（随机抽样 10 万条观测）：

下图展示了上市公司周边平均 NDVI 的数值分布：

下图按距离圈层展示了平均 NDVI 的分布（分办公地址 / 注册地址），可以清晰看到从公司所在地向外、植被指数随距离上升的「城市核心—郊区」梯度：

处理方法

原始上市公司地址数据来自 2000～2025年上市公司注册地址与办公地址（含经纬度、所处的省市区县及搬迁距离）.dta，栅格底图来自逐年合成的 年度tif/ 文件夹（包含 annual_mean 年度均值合成与 annual_max 年度最大值合成两套栅格）。处理流程分为三步，对应 code/ 文件夹中的代码。

第一步，使用 Stata 按年份拆分上市公司地址数据，便于后续逐年的并行处理：

*- 读取上市公司注册地址与办公地址数据（含经纬度）
use "2000～2025年上市公司注册地址与办公地址（含经纬度、所处的省市区县及搬迁距离）.dta", clear

*- 按年份拆分数据到「上市公司分年」文件夹
cap mkdir "上市公司分年"
split_by_var 年份, folderpath("上市公司分年") prefix("")

第二步，使用 R 语言（sf + terra）围绕每家公司的办公地址与注册地址，生成 10 个距离圈层的缓冲区分幅，并利用并行计算加速：

# 生成「年份 × 距离」组合
2001:2025 %>%
  crossing(c(1, 2, 3, 4, 5, 10, 15, 20, 25, 30)) %>%
  set_names("year", "dist") -> yddf

# 并行集群（每个节点加载 tidyverse 与 sf）
library(parallel)
makeCluster(5) -> cl
clusterExport(cl, "yddf")
clusterEvalQ(cl, ({
  library(tidyverse)
  library(sf)
}))

# 按年份读取地址数据，构造点并生成指定半径的缓冲区，写出为 rds
dir.create("bufferdf")
parLapply(cl, 1:nrow(yddf), function(x) {
  # 注册地址缓冲区
  haven::read_dta(paste0("上市公司分年/data_", yddf$year[x], ".dta")) %>%
    filter(!is.na(注册地址_经度)) %>%
    st_as_sf(coords = c("注册地址_经度", "注册地址_纬度"), crs = 4326) %>%
    st_buffer(dist = units::set_units(yddf$dist[x], "km")) %>%
    write_rds(paste0("bufferdf/", yddf$year[x], "_注册地址_", yddf$dist[x], ".rds"))

  # 办公地址缓冲区
  haven::read_dta(paste0("上市公司分年/data_", yddf$year[x], ".dta")) %>%
    filter(!is.na(办公地址_经度)) %>%
    st_as_sf(coords = c("办公地址_经度", "办公地址_纬度"), crs = 4326) %>%
    st_buffer(dist = units::set_units(yddf$dist[x], "km")) %>%
    write_rds(paste0("bufferdf/", yddf$year[x], "_办公地址_", yddf$dist[x], ".rds"))
}) -> res

第三步，使用 terra::extract() 从年度均值 / 最大值合成栅格中，分别提取每个缓冲区内的平均 NDVI 与最大 NDVI，并合并为面板数据写出：

# 逐缓冲区文件提取 NDVI（均值 / 最大值）
cl <- makeCluster(16)
clusterEvalQ(cl, ({
  library(tidyverse)
  library(sf)
  library(terra)
}))

parLapply(cl, 1:nrow(yddf), function(x) {
  readr::read_rds(yddf$value[x]) -> dfbuffer1

  # 年度均值合成 tif：提取缓冲区内的均值 NDVI
  terra::rast(paste0("年度tif/", yddf$year[x], "_annual_mean.tif")) -> rst_mean
  terra::extract(rst_mean, terra::vect(dfbuffer1),
                 fun = "mean", na.rm = TRUE, exact = F) %>%
    as_tibble() %>% set_names(c("ID", "平均NDVI")) -> ex_mean

  # 年度最大值合成 tif：提取缓冲区内的若干最大值 NDVI
  terra::rast(paste0("年度tif/", yddf$year[x], "_annual_max.tif")) -> rst_max
  terra::extract(rst_max, terra::vect(dfbuffer1),
                 fun = "max", na.rm = TRUE, exact = F) %>%
    as_tibble() %>% set_names(c("ID", "最大NDVI")) -> ex_max

  # 合并并写出单年单距离结果
  dfbuffer1 %>% st_drop_geometry() %>%
    mutate(ID = row_number(),
           地址类型 = yddf$地址类型[x],
           距离_km = yddf$距离_km[x]) %>%
    left_join(ex_mean, by = "ID") %>%
    left_join(ex_max, by = "ID") %>%
    select(-ID) %>%
    write_rds(paste0("res/", x, ".rds"))
})

stopCluster(cl)

# 合并所有结果并整理为面板数据
fs::dir_ls("res", recurse = TRUE, regexp = "[.]rds$") %>% map_dfr(read_rds) -> df
df %>%
  mutate(年份 = as.integer(年份)) %>%
  select(股票代码, 股票简称, 年份, 地址类型, 距离_km, 平均NDVI, 最大NDVI) %>%
  arrange(股票代码, 年份, 地址类型, 距离_km) -> df_out

haven::write_dta(df_out, "2001～2025 年上市公司周边的 NDVI 面板数据.dta",
                 label = "数据处理：微信公众号 RStata")

整个过程实际上就是提取公司周边圆圈范围内的栅格区域求平均或最大值，圆圈边缘按照比例计算进入圈内的部分。

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

数据引用格式

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

RStata 数据中心: 2001～2025 年上市公司周边的 NDVI 面板数据. 2026. https://tidyfriday.cn/rsdb2/

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

RStata Data Center: Panel Data on NDVI around Listed Companies, 2001–2025. 2026. https://tidyfriday.cn/rsdb2/

关联课程/数据推荐

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

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

R 语言计算 NDVI 时序数据变异系数栅格：https://rstata.duanshu.com/#/course/e9e25461e95145c88aa739f3dcd14368

中国各省市碳排放量是如何计算的？R 语言栅格数据转面板数据：https://rstata.duanshu.com/#/course/75de598dcfcc4ad0b9ff28bff27f6b83

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

点击这里跳转到 RStata 短书平台获取附件：2001～2025 年上市公司周边的 NDVI 面板数据

数据概览

图表展示

处理方法

数据引用格式

关联课程/数据推荐

相关文章

评论