library(tidyverse) library(sf)
read_sf("2019行政区划/县.shp") %>% rename(县 = NAME, 县代码 = PAC) %>% st_transform(4326) %>% dplyr::filter(!is.na(县) & 省 != "台湾省") %>% st_make_valid() -> county
st_intersects(st_make_valid(county), st_make_valid(county), sparse = F) -> touchmat touchmat %>% as_tibble() %>% mutate_all(as.numeric) %>% mutate(县代码 = county$县代码) %>% select(县代码, everything()) %>% set_names(c("县代码", county$县代码)) %>% gather(2:ncol(.), key = "key", value = "value") %>% dplyr::filter(县代码 != key & value == 1) -> pairdf
pairdf
pairdf %>% mutate(key = as.numeric(key)) %>% dplyr::filter(县代码 > key) %>% rename(相邻县 = key) %>% select(-value) -> pairdf
pairdf %>% slice(1) -> pair1
county %>% dplyr::filter(县代码 == pair1$县代码) -> maincounty county %>% dplyr::filter(县代码 == pair1$相邻县) -> touchcounty
st_intersection(st_cast(maincounty$geometry, "MULTILINESTRING"), st_cast(touchcounty$geometry, "MULTILINESTRING")) %>% st_collection_extract("LINESTRING") -> commonline
bind_rows(maincounty, touchcounty) %>% select(-contains("类型")) %>% st_centroid() -> twocounty
twocounty %>% mutate( 距离相邻县对共同边界的最小距离 = twocounty %>% st_distance(commonline) %>% units::set_units("km") ) %>% mutate(距离相邻县对共同边界的最小距离 = as.numeric(距离相邻县对共同边界的最小距离)) -> countycentroid
countycentroid %>% st_drop_geometry() %>% mutate(pair = paste0(县代码[1], "-", 县代码[2]))
ggplot() + geom_sf(data = maincounty, aes(fill = 市), alpha = 0.3) + geom_sf(data = touchcounty, aes(fill = 市), alpha = 0.5) + geom_sf(data = commonline, color = "#709ae1", size = 2) + geom_sf(data = countycentroid, aes(size = 距离相邻县对共同边界的最小距离), color = "gray20", shape = 16) + scale_fill_manual(values = c("#ff847c", "#99b898"), name = "") + scale_size_continuous(range = c(1, 4), name = "距离相邻县对共同边界\n的最小距离(km)") + labs(title = paste0(maincounty$县, "与", touchcounty$县, "质心与共同边界的最小距离"), subtitle = "数据计算&绘图:微信公众号 RStata", caption = "注:散点大小表示质心距离两个共同区界的最小距离,图中的蓝色粗线表示两个区的共同市界") + scale_x_continuous(limits = c(116.3, 116.52)) + theme(legend.position.inside = c(0.9, 0.5), legend.position = "inside")
|
评论