之前给大家分享过两个课程：
×
×
不过课程里面的南北方区域是手动绘制的，不是很方便，今天再给大家介绍一种新的方案。
首先加载所需的 R 包：
library(tidyverse)
library(sf)
library(lwgeom)
|
读取 2014 年工企地理位置数据：
haven::read_dta("2014年工企数据库地理位置数据.dta") %>%
mutate(id = row_number()) %>%
select(id, ends_with("度")) -> df
df
|
读取省级行政区划：
read_sf("2021行政区划/省.shp") -> prov
prov %>%
st_union() %>%
nngeo::st_remove_holes() -> cn
|
中国的经纬度范围：
st_bbox(cn) %>%
st_as_sfc() -> cnbbox
|
读取秦岭-淮河线与胡焕庸线
# 读取秦岭-淮河线
read_sf("秦岭淮河线/秦岭淮河线.shp") -> qh
# 读取胡焕庸线
read_sf("胡焕庸线/胡焕庸线.shp") -> hhy
|
我们再使用 leaflet 一并预览下上面的对象们：
library(leaflet)
library(mapview)
leaflet() %>%
addTiles("http://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineCommunity/MapServer/tile/{z}/{y}/{x}", attribution = "微信公众号 RStata") -> map
mapview(cnbbox, layer.name = "中国范围",
col.regions = "red", alpha.regions = 0.4) +
mapview(cn, layer.name = "中国") +
mapview(qh, color = "red", layer.name = "秦岭-淮河线") +
mapview(hhy, map = map, color = "blue", layer.name = "胡焕庸线")
|
如何使用胡焕庸线分东西
lwgeom::st_split() 函数可以线条把一个多边形分成两个区域，但是胡焕庸线的长度不够，所以我们得把它往两端延伸，这里我们可以选择 cnbbox 左下角和右上角的两个点。
# 胡焕庸线的两个端点
st_startpoint(hhy) -> hhytr
st_endpoint(hhy) -> hhybl
|
cnbbox 左下角和右上角的两个点：
cnbbox %>%
st_cast("POINT") %>%
.[1] -> blpoint
cnbbox %>%
st_cast("POINT") %>%
.[3] -> trpoint
|
把这四个点合成一条新的 linestring：
c(blpoint, hhybl, hhytr, trpoint) %>%
st_coordinates() %>%
st_linestring() %>%
st_sfc(crs = 4326) -> longhhyline
|
最终的效果：
longhhyline %>%
mapview() +
mapview(cnbbox, map = map)
|
这里还可以编写一个函数用于 point 合成 linestring：
st_point_to_linestring <- function(..., crs) {
c(...) %>%
st_coordinates() %>%
st_linestring() %>%
st_sfc(crs = crs)
}
st_point_to_linestring(blpoint, hhybl, hhytr, trpoint,
crs = 4326)
|
然后就可以进行分隔了：
cnbbox %>%
lwgeom::st_split(longhhyline) %>%
st_collection_extract("POLYGON") %>%
st_as_sf() %>%
mutate(xpos = st_coordinates(st_centroid(.))[,"X"]) %>%
mutate(position = ifelse(xpos == max(xpos), "east", "west")) -> west_eastdf
west_eastdf
|
上面代码中的 xpos 用于判断两个部分的东西情况。
mapview(west_eastdf, map = map, zcol = "position")
|
然后就可以对工企进行东西分类了。
首先把工企地理位置数据转换成 sf 对象：
df %>%
filter(!is.na(经度)) %>%
st_as_sf(coords = c("经度", "纬度"),
crs = 4326, remove = F) -> dfsf
dfsf %>%
st_intersection(west_eastdf) -> dfsf1
dfsf1
|
再计算工企距离胡焕庸线的距离：
dfsf1 %>%
st_distance(hhy) -> distmat
dfsf1 %>%
st_drop_geometry() %>%
mutate(dist = distmat[,1],
dist = units::set_units(dist, "km"),
dist = as.numeric(dist)) %>%
select(-xpos) -> dfsf1
|
如何使用秦岭-淮河线分南北
可以在 cnbbox 的左右框线取中点对秦岭淮河线进行延伸：
cnbbox %>%
st_cast("POINT") %>%
.[c(4, 5)] %>%
st_point_to_linestring(crs = 4326) %>%
st_centroid() -> leftpoint
cnbbox %>%
st_cast("POINT") %>%
.[c(2, 3)] %>%
st_point_to_linestring(crs = 4326) %>%
st_centroid() -> rightpoint
c(rightpoint, st_cast(qh$geometry, "POINT"), leftpoint) %>%
st_point_to_linestring(crs = 4326) -> longqhline
|
然后就可以进行分隔了：
cnbbox %>%
lwgeom::st_split(longqhline) %>%
st_collection_extract("POLYGON") %>%
st_as_sf() %>%
mutate(ypos = st_coordinates(st_centroid(.))[,"Y"]) %>%
mutate(position = ifelse(ypos == max(ypos), "north", "south")) -> north_southdf
|
绘图展示下：
mapview(north_southdf, map = map, zcol = "position")
|
然后就可以对工企进行南北分类了：
dfsf %>%
st_intersection(north_southdf) -> dfsf2
dfsf2
|
再计算工企距离秦岭-淮河线的距离：
dfsf %>%
st_distance(qh) -> distmat2
dfsf2 %>%
st_drop_geometry() %>%
mutate(dist = distmat2[,1],
dist = units::set_units(dist, "km"),
dist = as.numeric(dist)) %>%
select(-ypos) -> dfsf2
dfsf2
|
点击这里跳转到 RStata 短书平台获取附件：R 语言：如何根据经纬度判断所处的南北方、东西部
评论