最近有小伙伴问到了关于如何计算一些县到另外一些县质心最小距离的问题。今天我们一起来看下如何在 R 语言以及 Stata 中进行这样的计算。
首先我们使用 Stata 生成一些样本：
use county_db.dta, clear
keep if inlist(省, "江苏省", "安徽省", "浙江省", "山东省")
sample 10
keep 县 县代码
ren 县 试点县
export excel using 政策试点县.xlsx, replace first(variables)
use county_db.dta, clear
keep if inlist(省, "江苏省", "安徽省", "浙江省", "山东省")
keep 县 县代码
ren 县 样本县
export excel using 样本县.xlsx, replace first(variables)
|
我们的目标就是计算所有的样本县与政策试点县的最小距离。
使用 R 语言的 sf 包
R 语言中通常使用 sf 包进行矢量数据的地理计算。
首先加载相关 R 包，读取数据：
library(tidyverse)
library(sf)
read_sf("2020行政区划/县.shp") %>%
select(-contains("类型")) %>%
mutate(县代码 = as.numeric(县代码)) -> county
readxl::read_xlsx("政策试点县.xlsx") -> cty1
cty1
readxl::read_xlsx("样本县.xlsx") -> cty2
cty2
|
由于 cty1 和 cty2 中都没有地理信息，所以需要和 county 匹配下：
cty1 %>%
left_join(county, by = join_by("县代码")) %>%
st_sf() %>%
st_centroid() -> cty1
cty2 %>%
left_join(county, by = join_by("县代码")) %>%
st_sf() %>%
st_centroid() -> cty2
|
st_centroid() 函数的作用是计算每个区县的质心。
然后就可以使用 st_distance() 计算距离矩阵了：
cty2 %>%
st_distance(cty1) -> distmat
dim(distmat)
|
由于计算得到的距离矩阵中并没有县的信息，所以我们还得准备一个样本县的信息：
cty2 %>%
st_drop_geometry() %>%
select(样本县, 县代码) -> ctydf2
|
然后就可以求距离矩阵每一行的最小值，再和上面的样本县合并了：
apply(distmat, 1, min) %>%
as_tibble() %>%
bind_cols(ctydf2, .) %>%
mutate(value = value / 1000) %>%
rename(dist_km = value) -> resdf
resdf
resdf %>%
haven::write_dta("计算结果.dta")
|
这里计算的就是大圆距离。sf 包的计算规则是，如果计算的对象是 WGS84 坐标系，返回的是大圆距离，如果不是，返回的就是欧几里得距离。
使用 Stata 的 geodist 命令
Stata 的 geodist 命令计算的也是大圆距离。
首先读取两个数据：
import excel using "政策试点县.xlsx", clear first
save cty1, replace
import excel using "样本县.xlsx", clear first
foreach i of varlist _all {
cap ren `i' `i'2
}
save cty2, replace
|
注意，由于这两个数据的变量名有重复，所以为了下面的 corss 操作，我把 cty2 数据的变量名都加了个 2。
Stata 中可以使用 shp2dta 命令把 shp 格式的矢量数据转换成 dta 数据，使用 gencentroids(centroids) 选项可以计算出每个区域的质心。
*- 计算区县质心
shp2dta using 2020行政区划/县.shp, database(county_db) coordinates(county_coord) genid(ID) replace gencentroids(centroids)
|
然后得到的 county_db.dta 数据里面就有质心坐标 x_centroids（经度）和 y_centroids（纬度）了。
下面我们把 cty1 和 cty2 cross 起来，并且匹配上每个县的质心：
use cty2, clear
cross using cty1
merge m:1 县代码 using county_db
drop if _m == 2
drop _m
ren x_centroids lon1
ren y_centroids lat1
drop ID 省 - 县
drop 县代码
ren 县代码2 县代码
merge m:1 县代码 using county_db
drop if _m == 2
drop _m
ren x_centroids lon2
ren y_centroids lat2
drop ID 省 - 县
|
然后就可以使用 geodist 命令计算距离了：
geodist lat1 lon1 lat2 lon2, gen(dist_km2)
collapse (min) dist_km, by(县代码 样本县2)
ren 样本县2 样本县
|
R 和 Stata 的对比使用
也可以把 R 语言的结果也拷贝过来对比下：
merge 1:1 样本县 县代码 using 计算结果
ren dist_km R语言计算结果
ren dist_km2 Stata计算结果
drop _m
list in 1/10
|
结果差异很小～ 也可以画图对比下：
tw sc Stata计算结果 R语言计算结果 || lfit Stata计算结果 R语言计算结果, ///
xti("R 语言计算结果") yti("Stata 计算结果") ///
ti("如何计算区县质心的地理距离：R 语言和 Stata 的对比使用") ///
subti("数据处理&计算：微信公众号 RStata") ///
leg(off)
gr export pic1.png, width(4800) replace
|
点击这里跳转到 RStata 短书平台获取附件：如何计算区县质心的地理距离：R 语言和 Stata 的对比使用
评论