旧版本｜2007～2016 年税调企业周边的 NDVI 平均值与最大值

新版本在这里：https://rstata.duanshu.com/#/brief/course/f4fa9236fe724f7585e3e0e27a4467f8

之前给大家分享过 1982～2022 年 GIMMS NDVI 指数（归一化植被指数）各省市区县、乡镇面板数据：

1982～2022 年 GIMMS NDVI 指数（归一化植被指数）各省市区县、乡镇面板数据：https://rstata.duanshu.com/#/brief/course/fe5cc74d2a0d4dfcbb06d19fc99a30cf

最近有小伙伴表示想要计算历年税调企业周边一定范围内的 NDVI 数据。

税调企业的经纬度数据可以使用之前分享的：

2007~2016 年税调企业地理信息数据（含经纬度及其所处的省市区县）：https://rstata.duanshu.com/#/course/76d38022cd004b09b2aa09647936beb0

数据概览

为了方便大家使用，我计算了税调企业周边 1km、2km、3km、4km、5km 范围内的平均和最大年度 NDVI 数据，最终得到了附件中的：2007～2016 年税调企业周边的 NDVI 平均值与最大值.dta，数据概览：

包含如下变量：

年份、sdid、sdgroup、法人代码、企业名称、dist1km_max、dist1km_mean、dist2km_max、dist2km_mean、dist3km_max、dist3km_mean、dist4km_max、dist4km_mean、dist5km_max、dist5km_mean

下图展示了各年税调企业周边 NDVI 最大值的变化：

各年税调企业周边 NDVI 平均值的变化：

处理方法

首先我们可以使用类似这样的 R 语言代码构造税调企业周边 5km 的缓冲区：

df2 %>%
  st_buffer(dist = units::set_units(5, km)) -> df2_5km

然后再使用该缓冲区从 NDVI 栅格数据中提取该范围的 NDVI 均值或者最大值：

terra::extract(rst, vect(df2_1km), fun = y, na.rm = T, exact = T) %>%
    as_tibble() %>%
    mutate(year = y) %>%
    write_rds(paste0("计算结果/df2_1km_", y , ".rds"))

附件中也提供了详细的代码以供参考，对该数据计算感兴趣的小伙伴可以学习这个课程：

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

为了便于大家理解这个计算过程，我绘制了一幅这样的图：

方块表示的是 NDVI 栅格数据，圆圈表示的是税调企业周边 5km 范围。计算过程实际上就是提取圆圈范围内的区域求平均或者最大值。圆圈边缘的按照比例计算进入圈内的部分。

注意这里使用的栅格数据是年度生长季均值合成结果，计算结果中的均值是指区域的均值，最大值是指区域的最大值。

点击这里跳转到 RStata 短书平台获取附件：旧版本｜2007～2016 年税调企业周边的 NDVI 平均值与最大值