用 sf: st_intersection 分割线或多边形不起作用？

Question

我正在尝试完成沿乌拉尔山脉分裂俄罗斯这一看似简单的任务，但我一定是误解了 sf_intersection 中的一些基本内容。我什至无法分割线，更不用说创建两个新的多边形

rm(list=ls())
library(rworldmap)
library(tidyverse)
library(sf)

setwd("/Users/mr56267/Documents/UT_web/TLAH_Maps_2023/Mapping_textbook")
world_sp <- fortify(getMap())




world_sf <- world_sp %>% st_as_sf(coords = c("long", "lat"), crs = 4326, row.names="group") %>%
  group_by(group) %>% summarise(geometry = st_combine(geometry)) %>% st_cast("POLYGON")

ggplot() +
  geom_sf(data = world_sf)

world_robinson <- st_transform(world_sf, 
                               crs = '+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')



Russia.sf <- world_robinson$geometry[which(world_robinson$group=="Russia.1")]
### good to here
### now split Russia on Urals

https://raw.githubusercontent.com/sashatrubetskoy/asia_europe_border/master/asia_europe_border.geojson
Urals.sf <- st_read("https://raw.githubusercontent.com/sashatrubetskoy/asia_europe_border/master/asia_europe_border.geojson")
Urals.sf <- st_transform(Urals.sf, crs = st_crs(world_robinson$geometry[which(world_robinson$group=="Russia.1")]))
Russia_line.sf <- st_cast(Russia.sf,"MULTILINESTRING")

intersection_points <- st_intersection(Russia_line.sf,Urals.sf)
intersection_points
st_intersects(intersection_points,Russia_line.sf)

这里发生了什么？用 st_intersection 创建的点怎么可能不是交点呢？一切看起来都不错：

ggplot() + geom_sf(data=Russia_line.sf) + geom_sf(data=Urals.sf, color="red") + 
  geom_sf(data=intersection_points, color="green")

如何将俄罗斯分为东、西两部分？

Answer 1

按照评论中的建议使用lwgeom::st_split()，形式为lwgeom::st_split(POLYGON, LINESTRING). 需要记住的一件事是，您使用的边界非常详细，因此根据分辨率Russia.sf（例如边界的详细程度），它可能会产生不需要的伪影，特别是在里海附近（预计会有限制）以匹配陆地边界）。

giscoR查看使用具有合理良好分辨率的世界地图的完整示例（我进行了测试rworldmap，这就是出现工件的时候）：

library(giscoR)
library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE

# Using giscoR instead
world_sf <- gisco_get_countries(resolution = 10)


ggplot() +
  geom_sf(data = world_sf)

world_robinson <- st_transform(world_sf,
  crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
)

# Extract Russia and explode it in all the contiguous landmasses (POLYGONS)
Russia_exploded <- world_robinson %>%
  filter(ISO3_CODE == "RUS") %>%
  st_cast("POLYGON")
#> Warning in st_cast.sf(., "POLYGON"): repeating attributes for all
#> sub-geometries for which they may not be constant

# Extract mainland as the largest contiguous portion of Russia
# Equivalent to your world_robinson$group=="Russia.1"
Russia.sf <- Russia_exploded[which.max(st_area(Russia_exploded)), ] %>%
  st_union() %>%
  st_sf(cntry = "Russia")

ggplot(Russia.sf) +
  geom_sf()

### now split Russia on Urals
# https://raw.githubusercontent.com/sashatrubetskoy/asia_europe_border/master/asia_europe_border.geojson
download.file("https://raw.githubusercontent.com/sashatrubetskoy/asia_europe_border/master/asia_europe_border.geojson",
  "border.geojson",
  mode = "wb"
)

Urals.sf <- read_sf("border.geojson") %>%
  st_geometry() %>%
  st_combine() %>%
  st_line_merge() %>%
  st_transform(st_crs(world_robinson))


ggplot(Russia.sf) +
  geom_sf() +
  geom_sf(data = Urals.sf)

# Split with lwgeom
Russia_splitted <- lwgeom::st_split(Russia.sf, Urals.sf) %>%
  st_collection_extract("POLYGON")

Russia_splitted$gr <- factor(seq_len(nrow(Russia_splitted)))


ggplot(Russia_splitted) +
  geom_sf(aes(fill = gr))

# Done, now identify resulting polygons (visually) and assign continent
Russia_splitted_end <- Russia_splitted %>%
  mutate(Continent = ifelse(gr == "1", "Asia", "Europe")) %>%
  group_by(Continent) %>%
  summarise()

ggplot(Russia_splitted_end) +
  geom_sf(aes(fill = Continent))

# Same but faceted

ggplot(Russia_splitted_end) +
  geom_sf(aes(fill = Continent)) +
  facet_wrap(~Continent, ncol = 1)

^{创建于 2023-06-28，使用reprex v2.0.2}