之前给大家分享过历年工企周边的 PM2.5 浓度数据，今天我们一起来看下这种数据是如何使用 R 语言计算的。
历年工企业周边年平均 PM2.5 浓度面板数据：https://rstata.duanshu.com/#/brief/course/96e2daf6e4d94af9baf7f98908d6c535
逐年计算
首先我在附件中给大家准备了大概 1 万条工企数据的样本：
library(tidyverse)
library(sf)
library(terra)
haven::read_dta("工企样本.dta") -> df
df %>%
mutate(year = as.numeric(str_sub(as.character(gqid), 1, 4))) -> df
df
|
转换成 sf 对象：
df %>%
fill(经度, 纬度, .direction = "updown") %>%
st_as_sf(coords = c("经度", "纬度"), crs = 4326) -> dfsf
dfsf
|
首先以 2004 年为例讲解计算过程：
dfsf %>%
filter(year == 2004) -> dfsf2004
dfsf2004
|
生成坐标点周边 5km 的缓冲区（圆形区域）：
dfsf2004 %>%
st_buffer(dist = units::set_units(5, km)) -> dfbuffer_5km
|
这个区域数据是这样的：
library(leaflet)
url <- "http://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineStreetWarm/MapServer/tile/{z}/{y}/{x}"
leaflet() %>%
addTiles(url, attribution = "微信公众号 RStata") -> map
dfbuffer_5km %>%
filter(str_detect(企业名称, "北京")) %>%
mapview::mapview(map = map)
|
读取 2004 年的栅格数据：
rast("yeartif/cn2004.tif") -> rst
rst
#> class : SpatRaster
#> dimensions : 4724, 6159, 1 (nrow, ncol, nlyr)
#> resolution : 0.01, 0.01 (x, y)
#> extent : 73.5, 135.09, 6.320002, 53.56 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source : cn2004.tif
#> name : cn2004
#> min value : 1
#> max value : 108
|
然后就可以提取各工企 5km 范围内的 PM2.5 浓度均值了：
terra::extract(rst, vect(dfbuffer_5km),
fun = "mean", na.rm = T, exact = T) %>%
as_tibble() -> res5km
|
再和原数据合并：
dfbuffer_5km %>%
st_drop_geometry() %>%
mutate(ID = row_number()) %>%
left_join(res5km) %>%
select(-ID) %>%
rename(PM25_5km = cn2004) -> df5km
df5km
|
对于多年数据，循环各年然后合并即可：
dfsf %>%
st_buffer(dist = units::set_units(5, km)) -> dfsfbuffer_5km
dfsfbuffer_5km
lapply(1998:2014, function(x){
print(x)
dfsfbuffer_5km %>%
filter(year == x) -> buffertemp
rast(paste0("yeartif/cn", x, ".tif")) -> rsttemp
terra::extract(rsttemp, vect(buffertemp),
fun = "mean", na.rm = T, exact = T) %>%
as_tibble() -> restemp
buffertemp %>%
st_drop_geometry() %>%
mutate(ID = row_number()) %>%
left_join(restemp) %>%
select(-ID) %>%
set_names("group", "gqid", "企业名称",
"year", "PM25_5km")
}) %>%
bind_rows() -> dfall
dfall
|
这样就得到了所有企业各年周边 5km 范围内的平均 PM2.5 浓度。
如果想要多个距离的，可以再嵌套一个循环，例如计算 5 和 10km 的：
lapply(c(5, 10), function(y){
print(y)
dfsf %>%
st_buffer(dist = units::set_units(y, km)) -> dfsfbuffer_temp
lapply(1998:2014, function(x){
dfsfbuffer_temp %>%
filter(year == x) -> buffertemp
rast(paste0("yeartif/cn", x, ".tif")) -> rsttemp
terra::extract(rsttemp, vect(buffertemp),
fun = "mean", na.rm = T, exact = T) %>%
as_tibble() -> restemp
suppressMessages({
buffertemp %>%
st_drop_geometry() %>%
mutate(ID = row_number()) %>%
left_join(restemp) %>%
select(-ID) %>%
set_names("group", "gqid", "企业名称",
"year", paste0("PM25_", y, "km"))
})
}) %>%
bind_rows()
}) -> resall
resall %>%
reduce(left_join) -> resalldf
resalldf
resalldf %>%
haven::write_dta("1998～2014年部分工企周边的平均PM2.5浓度.dta")
|
绘图：
haven::read_dta("1998～2014年部分工企周边的平均PM2.5浓度.dta") -> resalldf
library(gghighlight)
ggplot(resalldf, aes(x = year, y = PM25_5km,
color = factor(group))) +
geom_line() +
gghighlight(group %in% sample(unique(resalldf$group),
size = 8),
label_params = list(family = cnfont,
label = "企业名称")) +
scale_color_brewer(palette = "Set1") +
ggthemes::theme_solarized_2(base_family = cnfont) +
theme(plot.background = element_rect(fill = "#fff5e3"),
plot.margin = grid::unit(rep(0.8, 4), "cm"),
axis.title.x = element_blank()) +
theme(legend.position = "none") +
scale_x_continuous(breaks = 1998:2014) +
scale_y_continuous(breaks = seq(0, 150, by = 20)) +
labs(y = "年均 PM2.5 质量浓度(µg/m3)",
title = "1998~2014年部分工企周边 5km 范围内的平均 PM2.5 浓度",
subtitle = "数据计算：微信公众号 RStata",
caption = "数据来源：工企数据库，经纬度根据高德地图地理编码接口解析，PM2.5 栅格数据来源于华盛顿大学圣路易斯分校") -> p
|
这个过程的计算原理可以通过下面的图理解：
library(raster)
read_sf("北京长条范围.geojson") -> fanwei
raster("yeartif/cn2004.tif") -> rst
rst %>%
raster::mask(fanwei) %>%
raster::crop(fanwei) %>%
rasterToPolygons() -> rstpoly
dfbuffer_5km %>%
st_intersection(fanwei) -> dftemp
rstpoly %>%
st_as_sf() %>%
st_intersection(fanwei) -> rstpolytemp
library(mapview)
mapview(rstpolytemp,
layer.name = "NDVI") +
mapview(dftemp,
layer.name = "5 km", map = map)
|
再回想前面的提取代码为：
也就是加总平均每个圆圈内的值，如果使用了 exact = T 参数，表示边缘的方块根据进入内部的比例计算，如果 exact = F，表示边缘的方块不作考虑，只计算内部的。
直接计算多年的
另外也可以使用这样的方法计算：
fs::dir_ls("yeartif")
rast(fs::dir_ls("yeartif")) -> rst
terra::extract(rst, vect(dfsfbuffer_5km),
fun = "mean", na.rm = T, exact = T) %>%
as_tibble() -> res5km
res5km
|
再合并原数据、去除不同年份的结果即可：
dfsfbuffer_5km %>%
st_drop_geometry() %>%
mutate(ID = row_number()) %>%
left_join(res5km) %>%
select(-ID) %>%
gather(contains("cn"), key = "key", value = "PM25_5km") %>%
mutate(key = str_remove_all(key, "cn"),
key = as.numeric(key)) %>%
filter(year == key) %>%
slice(-key)
|
有时候我们也会遇到不同年份的栅格数据范围不一致、分辨率不一致的问题（会无法直接使用 rast 读取），这个时候可以先把各年的栅格数据重采样成一致的，然后再使用 rast() 读取：
fs::dir_ls("yeartif") -> fls
rast(fls[1]) -> rst1
lapply(fls, function(x){
rast(x) %>%
terra::resample(rst1)
}) %>%
rast() -> rstnew
names(rstnew) <- 1998:2014
rstnew
|
然后后面的操作就都一样了。
点击这里跳转到 RStata 短书平台获取附件：使用 R 语言计算工企周边一定范围内的平均 PM2.5 浓度
评论