之前给大家分享过很多污染物数据，例如 SO2、粉尘等指标：
- 1980 年 1 月～2025 年 12 月各省市区县地表 SO2 质量浓度：https://rstata.duanshu.com/#/course/72351a45d60c4c63af3de00db732221c
- 1980 年 1 月～2025 年 12 月各省市区县地表粉尘质量浓度：https://rstata.duanshu.com/#/brief/course/6f1ed61d557e4d9aa388a8154320b801
- …（更多污染物数据可以在平台搜索）
最近有小伙伴表示想要计算历年税调企业周边的卫星反演污染物数据，考虑到计算量比较大，我就直接帮忙处理好了。
税调企业的经纬度数据可以使用之前分享的：
2007~2020 年税调企业地理信息数据（含经纬度及其所处的省市区县）：https://rstata.duanshu.com/#/course/39f5d95ad58841628fa4d935371bf3fd
数据概览
为了方便大家使用，我计算了税调企业周边 1km、2km、3km、4km、5km 范围内各种污染物指标的年度平均数据（因为计算量非常大，所以这次没有计算 10km、15km 和 20km 的）。时间范围为 2007～2020 年，经过计算之后一共得到了 819.9 万条结果：
这里我提供的数据格式为 dta 格式，包含的变量如下：
年份、sdgroup、企业名称、法人代码、sdid、PM1_1km、PM1_2km、PM1_3km、PM1_4km、PM1_5km、PM10_1km、PM10_2km、PM10_3km、PM10_4km、PM10_5km、PM2_5_1km、PM2_5_2km、PM2_5_3km、PM2_5_4km、PM2_5_5km、SO2_1km、SO2_2km、SO2_3km、SO2_4km、SO2_5km、SO4_1km、SO4_2km、SO4_3km、SO4_4km、SO4_5km、二甲硫醚_1km、二甲硫醚_2km、二甲硫醚_3km、二甲硫醚_4km、二甲硫醚_5km、粉尘_1km、粉尘_2km、粉尘_3km、粉尘_4km、粉尘_5km、黑碳_1km、黑碳_2km、黑碳_3km、黑碳_4km、黑碳_5km
绘图展示
为了更直观地感受这份数据，我还绘制了两幅图进行展示。下图展示了历年税调企业周边 1～5km 范围内年平均 PM10 浓度：
SO2 浓度：
处理方法
处理方法类似之前的课程「使用 R 语言计算工企周边一定范围内的平均 PM2.5 浓度」，感兴趣的小伙伴可以学习：
使用 R 语言计算工企周边一定范围内的平均 PM2.5 浓度：https://rstata.duanshu.com/#/course/016330182fb64ecd9fab8ada0d816130
具体处理方法如下：
- 首先分别生成税调企业地址周边 1km、2km、3km、4km、5km 的缓冲区。
- 读取各种污染物指标的栅格数据。
- 使用上面的缓冲区数据和读取到的栅格数据进行地理计算，就可以得到每个缓冲区内的各种污染物指标的均值了。为了便于大家理解这个计算过程，我绘制了一幅下面的图，方块表示的是污染物栅格数据，圆圈表示的是税调企业周边 5km 范围。计算过程实际上就是提取圆圈范围内的区域求平均值。圆圈边缘的按照比例计算进入圈内的部分。
按照这个思路循环各个年份的即可（不过运算量非常大，建议大家不要尝试，非常耗时）。
处理之后再整理即可得到分享给大家的税调企业周边的卫星反演年度污染物数据.dta了。
附件中也提供了匹配使用的代码供大家参考。
点击这里跳转到 RStata 短书平台获取附件：2007～2020 年税调企业周边的卫星反演年度污染物数据
评论