历年工企业周边年平均 PM2.5 浓度面板数据
前几天给大家分享了工企周边 1～5km 范围内的碳排放量数据：1998～2014 年工企周边 1～5km 范围内碳排放总量面板数据
很多小伙伴提出是不是可以使用类似的方法计算工企周边的 PM2.5 浓度数据。
之前给大家分享过不少关于 PM2.5 的数据：
日度数据：1980年1月1日～2022年6月30日中国各省市区县地表PM2.5质量浓度：https://rstata.duanshu.com/#/brief/course/ed262631255646ab856f6b2c970f99f8 卫星反演数据：1980 年 1 月～2022 年 3 月中国各省市区县地表 PM2.5 质量浓度面板数据（微克/立方米）：https://rstata.duanshu.com/#/brief/course/d96ee17d0fbc4fd78d11d41c002e0fca 2014 年至今中国各省市区县分年、分月、逐日 PM2.5 浓度面板数据：https://rstata.duanshu.com/#/brief/course/7a824c1de4004fcb9534a8c1fd902164 1998年1月～2019年12月中国各省、各城市、各区县 PM2.5 浓度月度面板数据 & 栅格数据（华盛顿大学圣路易斯分校）：https://rstata.duanshu.com/#/brief/course/3d48960e1cd240d2baf8b33135980500 1998～2019 年中国各省、各城市、各区县 PM2.5 浓度面板数据 & 栅格数据（华盛顿大学圣路易斯分校）：https://rstata.duanshu.com/#/brief/course/6a433349124e42d9b3cfb8d59a1df1e9 1998～2016 年中国各省、各城市、各区县 PM2.5 浓度面板数据 & 栅格数据（哥伦比亚大学）：https://rstata.duanshu.com/#/brief/course/69b7591a9f814d59b556c769c22bb5a7
关于这些数据是如何处理的，感兴趣的小伙伴可以参考这个课程：
中国各省市碳排放量是如何计算的？R 语言栅格数据转面板数据：https://rstata.duanshu.com/#/brief/course/75de598dcfcc4ad0b9ff28bff27f6b83
由于工企数据是 1998～2014 年间的，所以 1、2、4、5、6 都可以用，其中 1 和 2 是卫星反演数据，4、5 和 6 是结合卫星反演数据并使用地面观测站数据进行校正。所以我觉得用下面三个可能会更好，所以我这次先选择了第 5 份，华盛顿大学圣路易斯分校来源的 PM2.5 数据。
这份数据里面的栅格数据是分辨率 0.01˚x0.01˚ 的，大概是 1 km，下图展示了 2019 年的 PM2.5 浓度分布：
使用工企的经纬度坐标构建 1km～5km 的缓冲区就可以裁剪栅格数据计算各个工企历年周边的平均 PM2.5 浓度变化了。
工企数据我使用的是之前分享的这个：「1998～2014 年工业企业数据库地理位置数据（含经纬度、所处省市区县、南北方属性以及距离秦岭淮河线的距离）」（https://rstata.duanshu.com/#/brief/course/1c0d065e272c46b49cc528a445c2d3d9），其中经纬度是使用高德地图地理编码接口进行解析的。
首先分别生成了工企地址周边 1km、2km、3km、4km、5km 的缓冲区（因为计算量非常大，所以这次没有计算 10km、15km 和 20km 的），例如地址 5km 的缓冲区是这样生成的：
library(tidyverse)
library(sf)
haven::read_dta("1998～2014年工企地理位置面板数据.dta") %>%
filter(年份 == 2014) %>%
select(gqid, contains("度")) -> df1
df1 %>%
st_as_sf(coords = c("经度", "纬度"), crs = 4326) -> df2
df2 %>%
st_buffer(dist = units::set_units(5, km)) -> df2_5km
这样会得到一个以每个工企为中心、半径为 5km 的圆形区域。
然后读取碳排放量栅格数据，例如 2014 年的：
library(terra)
rast("yeartif/cn2014.tif") -> rst
然后使用上面的缓冲区数据和栅格数据进行地理计算就可以得到每个缓冲区内的碳排放总量了！例如 5km 的：
terra::extract(rst, vect(df2_5km), fun = "sum", na.rm = T, exact = T) %>%
as_tibble()
按照这个思路循环各个年份的即可（不过运算量非常大，建议大家不要尝试，非常耗时）。
处理之后再整理即可得到分享给大家的 历年工企周边的平均PM2.5浓度 数据了，数据中的浓度单位是“微克/立方米”：
例如北京宏华电器有限公司周边 1～5km 范围内年平均 PM2.5 浓度变化如下：
1～5km 的差异不大，大家可以根据自己的需要选择使用。
获取数据
是不是感觉很硬核！欢迎报名 RStata 培训班获取全部课程和以会员价获取数据资料（10元/份）详情可阅读这篇推文：数据处理、图表绘制、效率分析与计量经济学如何学习～
之前更新的课程和数据资料可点击阅读原文进入 RStata 学院查看（从首页的会员卡专区即可查看和购买会员卡）。
更多关于 RStata 培训班的信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/96e2daf6e4d94af9baf7f98908d6c535