历年工企周边 1 ～ 5 km范围内平均夜间灯光亮度

前几天给大家分享了工企周边 1～5km 范围内的碳排放量数据：

很多小伙伴提出是不是可以使用类似的方法计算工企周边的平均夜间灯光亮度数据。

之前给大家分享过不少关于夜间灯光的数据：

关于这些数据是如何处理的，感兴趣的小伙伴可以参考这个课程：

由于工企数据是 1998～2014 年间的，所以 2、3 都可以用，这次我先选择第二个进行计算。

这份数据里面的栅格数据是分辨率 1kmx1km 的，坐标系是 WGS_1984_Albers，下图展示了 2006 年的夜间灯光亮度分布：

使用工企的经纬度坐标构建 1km～5km 的缓冲区（注意也需要转换成 WGS_1984_Albers 坐标系）就可以裁剪栅格数据计算各个工企历年周边的平均夜间灯光亮度变化了。

工企数据我使用的是之前分享的这个：

其中经纬度是使用高德地图地理编码接口进行解析的。

首先分别生成了工企地址周边 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)) %>%
  st_transform("+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") -> df2_5km

这样会得到一个以每个工企为中心、半径为 5km 的圆形区域。

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

library(terra)
rast("yeartif/DMSP2014.tif") -> rst

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

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

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

处理之后再整理即可得到分享给大家的 工企周边的平均夜间灯光亮度数据.dta 数据了，：

例如常州电梯厂有限公司周边 1～5km 范围内年平均夜间灯光亮度变化如下：

1～5km 的差异不大，大家可以根据自己的需要选择使用。

点击这里跳转到 RStata 短书平台获取附件：历年工企周边 1 ～ 5 km范围内平均夜间灯光亮度

相关文章

评论