2007～2016 年税调周边 1～5km 范围内的碳排放总量（CGER来源）

前不久给大家分享过一份 2007～2016 年税调周边 1～5km 范围内的碳排放总量，其中碳排放栅格数据来源于 EDGAR（Emissions Database for Global Atmospheric Research）v8.0 版本：

2007～2016 年税调周边 1～5km 范围内的碳排放总量：https://rstata.duanshu.com/#/course/0fade2f56be342ea81092c7db266b603

今天再给大家分享一份碳排放数据来源于 CGER（Center for Global Environmental Research）的历年税调周边一定范围内碳排放总量的计算结果。

数据概览

同样地，我计算了税调企业周边 1km、2km、3km、4km、5km 范围内的碳排放总量，最终得到了附件中的：

• 2007～2016年税调周边1～5km范围内的碳排放总量（单位吨，CGER来源）.dta

经过计算之后一共得到了 622.7 万条结果，包含如下变量：

• 年份
• sdid
• sdgroup
• 法人代码
• 企业名称
• CO2_dist1km
• CO2_dist2km
• CO2_dist3km
• CO2_dist4km
• CO2_dist5km

下图展示了历年税调企业周边 1～5km 范围的的碳排放量（吨）：

处理方法

首先，分别生成税调企业注册地址周边 1km、2km、3km、4km、5km 的缓冲区，例如注册地址 5km 的缓冲区是这样生成的：

library(tidyverse)
library(sf)
haven::read_dta("2007~2016年税调经纬度及其所处的省市区县数据.dta") %>%
  filter(年份 == 2016) %>%
  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 的圆形区域。

然后读取碳排放量栅格数据，例如 2016 年的：

library(terra)
rast("cntif/2016.tif") -> rst

然后使用上面的缓冲区数据和栅格数据进行地理计算就可以得到每个缓冲区内的碳排放总量了！例如 5km 的：

terra::extract(rst, vect(df2_5km), fun = "sum", na.rm = T, exact = T) %>%
  as_tibble()

按照这个思路循环各个年份的即可（不过运算量非常大，建议大家不要尝试，非常耗时）。

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

方格子就是栅格数据，圆圈是税调企业周围 5km 范围的缓冲区。例如某个缓冲区和四个格子交叉，那么就是根据交叉部分的比例计算其范围内的碳排放总量。

附件中也提供了详细代码以供参考～

点击这里跳转到 RStata 短书平台获取附件：旧版本｜2007～2016 年税调周边 1～5km 范围内的碳排放总量（CGER来源）