在「城市空间结构与劳动者工资收入（刘修岩等，2019）」一文中，作者提出了三个描述城市区域在地面投影的平面几何形状特征的指标：
- road：城市中所有夜间灯光栅格之间地表直线距离的均值除以城市面积的平方根。反映城市内各地点之间的平均相互距离，代表绝大多数居民的出行成本。
- max：城市内所有夜间灯光栅格间距中的最大值除以城市面积的平方根。表征城市可能出现的最远距离，突出”形态不规则”或”多中心”特征。
- center：城市内每个夜间灯光栅格离城市中心点地表直线距离的均值除以城市面积的平方根。衡量城市中心的通达性，值越大说明大多数居民到市中心距离越远。
这三个指标从不同角度反映城市内部空间距离：road 关注平均距离，max 关注极端距离，center 关注中心可达性。指标值越高，城市形态越”劣质”。研究表明，城市形态通过影响居民出行成本和企业运输效率，对城市经济绩效和劳动者工资收入产生显著影响。
灯光阈值选择依据
夜间灯光数据记录了城市、城镇和乡村中相对稳定的灯光，年度夜光数据分辨率为 1kmx1km。作者选择亮度阈值40来界定城市化区域，基于以下考虑：
- 国际比较基准：Harari(2015) 对印度城市化区域的判定标准为亮度大于 35。考虑到中国的经济发展水平和能源消耗水平高于印度，我们适当提高阈值至40。
- 实证稳健性检验：使用 35 作为阈值重新计算城市形态指标后，研究结论保持一致，说明结果不依赖于特定阈值选择。
- 数据质量保证：阈值 40 能有效过滤乡村和小城镇的零星灯光，聚焦于真正的城市化区域。
本课程中将以 40 为例展示这几个指标的 R 语言计算过程。
数据处理流程与 R 代码详解
数据准备与重采样
夜间灯光数据使用的是平台上分享的这个：https://rstata.duanshu.com/#/brief/course/2c3880893f5947a5bf764befd3994c74，原始栅格数据可以从 Harvard dataverse 下载：https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/GIYGJU
这里我仅仅选择 2022～2024 年的数据演示：
library(tidyverse)
library(sf)
library(terra)
|
raw-data 里面存放的是 2022～2024 年的夜间灯光栅格数据，不过存在各年数据范围、分辨率不一致的问题，不方便进行批量运算，因此需要先通过重采样对齐：
fs::dir_ls("raw-data") -> ls
rast(ls[length(ls)]) -> rsttemplate
dir.create("newrst")
lapply(1:length(ls), function(x){
print(x)
rast(ls[x]) %>%
terra::resample(rsttemplate) %>%
writeRaster(paste0("newrst/", str_extract(ls[x], "\\d{4}"), ".tif"))
}) -> res
|
rast(fs::dir_ls("newrst")) -> rst
rst
|
按照作者的定义：城市区域为市辖区灯光亮度大于 40 的：
read_sf("2021行政区划/县.shp") %>%
st_make_valid() -> county
county %>%
filter(县类型 == "市辖区") %>%
group_by(省, 省代码, 市, 市代码) %>%
summarise() %>%
ungroup() -> city
city
library(terra)
fs::dir_ls("newrst") %>%
rast() -> rst
rst
city %>%
st_transform(crs(rst)) -> city
city
|
然后就可以分区域计算每个城市中夜光亮度大于等于 40 的面积了：
cellSize(rst, unit = "km") -> arearst
arearst[rst < 40] <- NA
names(arearst) <- names(rst)
plot(arearst)
|
可以看出在当前坐标系下各个像元的面积基本是一样的，不过我们还是直接加总每个城市区域里面的栅格面积，而不是仅仅统计数量：
terra::extract(arearst, vect(city),
fun = "sum", na.rm = T) %>%
as_tibble() -> res
bind_cols(city, res) %>%
select(-ID) %>%
st_drop_geometry() %>%
gather(contains("DMSP"), key = "year", value = "area") %>%
mutate(year = str_remove_all(year, "DMSP"),
year = as.numeric(year)) -> areadf
areadf
|
road、max、center 指标
然后就可以计算各年各城市的三个指标了：
mean_distance <- function(r, na.rm = T) {
dist_matrix <- terra::distance(as.points(r),as.points(r))
mean(dist_matrix[lower.tri(dist_matrix)], na.rm = T)
}
max_distance <- function(r, na.rm = T) {
dist_matrix <- terra::distance(as.points(r),as.points(r))
max(dist_matrix[lower.tri(dist_matrix)], na.rm = T)
}
rst[rst < 40] <- NA
rst
fc <- city[1,]
fc
rst %>%
crop(vect(fc)) %>%
mask(vect(fc)) -> frst
lapply(1:nlyr(frst), function(x){
mean_distance(frst[[x]], na.rm = T)
}) -> tempres
unlist(tempres)
frst[[1]] %>%
as.polygons() %>%
st_as_sf() %>%
st_union() %>%
st_cast("POLYGON") %>%
st_sf() %>%
mutate(area = st_area(.)) %>%
arrange(desc(area)) %>%
slice(1) %>%
st_centroid() -> centroids
centroids
plot(frst[[1]])
plot(centroids, add = T, col = "red")
|
然后就可以计算所有年份、所有城市的：
dir.create("res")
lapply(1:nrow(city), function(y){
fc <- city[y,]
rst %>%
crop(vect(fc)) %>%
mask(vect(fc)) -> frst
lapply(1:nlyr(frst), function(x){
if(!is.na(mean(frst[[x]][], na.rm = T))) {
frst[[x]] %>%
as.polygons() %>%
st_as_sf() %>%
st_union() %>%
st_cast("POLYGON") %>%
st_sf() %>%
mutate(area = st_area(.)) %>%
arrange(desc(area)) %>%
slice(1) %>%
st_centroid() -> centroids
frst[[x]] %>%
as.points() %>%
distance(vect(centroids)) -> distmat
fc %>%
st_drop_geometry() %>%
mutate(mean = mean_distance(frst[[x]], na.rm = T),
max = max_distance(frst[[x]], na.rm = T),
center = mean(distmat, na.rm = T) / 1000,
year = as.numeric(str_extract(names(frst)[x], "\\d{4}"))) %>%
write_rds(paste0("res/", y, "_", x, ".rds"))
}
}) -> tempres1
}) -> tempres2
|
合并计算结果：
fs::dir_ls("res") %>%
lapply(readr::read_rds) %>%
bind_rows() -> df2
df2
|
再合并两部分计算结果：
df2 %>%
full_join(areadf) %>%
mutate(road = mean / (1000 * sqrt(area)),
max = max / (1000 * sqrt(area)),
center = center / (sqrt(area))) %>%
select(-mean) %>%
select(省, 省代码, 市, 市代码, year, everything()) %>%
arrange(省, 省代码, 市, 市代码, year) %>%
rename(年份 = year) %>%
haven::write_dta("2022～2024年各城市形态指标面板数据.dta", label = "数据处理：微信公众号 RStata")
|
这里需要注意 mean 和 max 本来的单位是 m，所以除以 1000，center 我计算的单位就是 km 就不用再除以 1000 了。
点击这里跳转到 RStata 短书平台获取附件：R 语言：使用夜间灯光数据构造城市形态指标
评论