关于工企数据库和边界的距离，之前推出过不少很有意思的数据：
×
×
×
×
最近有个小伙伴表示自己需要工企距离省界、市界和县界的距离，于是就帮他算了下，这个计算过程挺简单，类似这个课程：
×
工企地理位置数据
首先可以从平台上下载到工企业的地理位置数据：
×
这里使用的是 高德地图结果合成面板 里面的那个：
首先合并：
use "1998.dta", clear
forval i = 1999/2014 {
di `i'
append using "`i'"
}
save "工企面板数据", replace
keep gqid 县代码 - 市 经度 纬度
save gqsim, replace
|
然后下面我们使用 R 语言计算距离：
library(tidyverse)
library(sf)
haven::read_dta("gqsim.dta") -> gqsim
gqsim %>%
dplyr::filter(!is.na(经度)) -> gqsim
gqsim %>%
st_as_sf(coords = c("经度", "纬度"), crs = 4326) -> gqsimsf
read_sf("2019行政区划/省.shp") %>%
st_transform(4326) %>%
dplyr::filter(省 != "中朝共有") %>%
st_cast("MULTILINESTRING") -> provline
read_sf("2019行政区划/市.shp") %>%
st_transform(4326) %>%
st_cast("MULTILINESTRING") -> cityline
read_sf("2019行政区划/县.shp") %>%
st_transform(4326) %>%
st_cast("MULTILINESTRING") -> countyline
|
然后就可以计算距离了，由于每个区域的点距离边界的最小距离肯定是该区域的边界，所以我们可以一个一个区域的计算，例如计算安徽省的：
gqsimsf %>%
dplyr::filter(省代码 == 340000) -> temp
provline %>%
dplyr::filter(省代码 == 340000) -> tempprov
temp %>%
st_distance(tempprov) -> dist
temp$dist <- dist[,1]
temp
|
那么我们只要循环所有的省份就可以得到全部的数据了，使用并行运算会更有效。
首先我们把需要的数据保存成 rds 文件：
gqsimsf %>%
write_rds("gqsimsf.rds")
provline %>%
write_rds("provline.rds")
cityline %>%
write_rds("cityline.rds")
countyline %>%
write_rds("countyline.rds")
|
然后运行 usethis::edit_r_profile() 打开 Profile 文件并把下面的代码复制进去保存后关闭：
library(tidyverse)
library(sf)
setwd("/Users/ac/Desktop/工企距离省市区县边界的距离/")
read_rds("gqsimsf.rds") -> gqsimsf
read_rds("provline.rds") -> provline
read_rds("cityline.rds") -> cityline
read_rds("countyline.rds") -> countyline
|
然后就可以根据自己的电脑性能创建多个集群：
library(parallel)
makeCluster(8) -> cl
|
并行处理省界：
parLapply(cl, provline$省代码, function(x){
gqsimsf %>%
dplyr::filter(省代码 == x) -> temp
provline %>%
dplyr::filter(省代码 == x) -> tempprov
temp %>%
st_distance(tempprov) -> dist
temp$dist <- dist[,1]
return(temp)
}) %>%
bind_rows() -> provdf
provdf %>%
write_rds("provdf.rds")
|
类似的方法处理市界和县界：
parLapply(cl, cityline$市代码, function(x){
gqsimsf %>%
dplyr::filter(市代码 == x) -> temp
cityline %>%
dplyr::filter(市代码 == x) -> tempcity
temp %>%
st_distance(tempcity) -> dist
temp$dist <- dist[,1]
return(temp)
}) %>%
bind_rows() -> citydf
citydf %>%
write_rds("citydf.rds")
parLapply(cl, countyline$PAC, function(x){
gqsimsf %>%
dplyr::filter(县代码 == x) -> temp
countyline %>%
dplyr::filter(PAC == x) -> tempcounty
temp %>%
st_distance(tempcounty) -> dist
temp$dist <- dist[,1]
return(temp)
}) %>%
bind_rows() -> countydf
countydf %>%
write_rds("countydf.rds")
|
然后合并三个文件即可：
provdf %>%
st_drop_geometry() %>%
select(gqid, dist) %>%
mutate(dist = as.numeric(dist),
dist = dist / 1000) %>%
rename(与省界的距离 = dist) %>%
left_join(
citydf %>%
st_drop_geometry() %>%
select(gqid, dist) %>%
mutate(dist = as.numeric(dist),
dist = dist / 1000) %>%
rename(与市界的距离 = dist)
) %>%
left_join(
countydf %>%
st_drop_geometry() %>%
select(gqid, dist) %>%
mutate(dist = as.numeric(dist),
dist = dist / 1000) %>%
rename(与县界的距离 = dist)
) -> df
df %>%
haven::write_dta("df.dta")
df
|
然后回到 Stata 里面再和原数据合并：
use 工企面板数据, clear
merge 1:1 gqid using df.dta
gsort gqid
drop _m
label var 与省界的距离 "与省界的最小距离"
label var 与市界的距离 "与市界的最小距离"
label var 与县界的距离 "与县界的最小距离"
label data "整理：微信公众号 RStata"
save "1998～2014年工企距离最近省界、市界、县界的距离", replace
|
由于这个文件高达 13.55GB，大家的电脑可能无法读取，所以我拆分下：
forval i = 1998/2014 {
use "1998～2014年工企距离最近省界、市界、县界的距离", clear
keep if 年份 == `i'
local N = `=_N'
foreach j of varlist _all {
qui count if missing(`j')
if r(N) == `N' {
drop `j'
}
}
save "data/`i'年工企距离最近省界、市界、县界的距离", replace
}
|
大家使用的时候可以逐年处理仅仅保存自己需要的变量再合并，这样合并后的文件就会小很多。
最后我绘制了一幅北京市的地图展示计算结果：
散点的颜色和大小表示该企业距离北京市边界的距离。
关于路网地图的绘制可以从平台上找到相关课程学习：
×
×
点击这里跳转到 RStata 短书平台获取附件：1998～2014 年工企与省份边界、城市边界及区县边界的最小距离（单位 km）
评论