问题

I'm trying to do the seemingly simple task of split Russia along the Ural mountains, but I must be misunderstanding something basic in sf_intersection. I can't even split the line, much less create two new polygons.

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)

What's going on here? How can the points created with st_intersection not be points of intersection? Everything looks good:

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

How do I split Russia into east and west?

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)

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

答案1

得分: 3

使用建议中提到的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

# 使用 giscoR
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"
)

# 提取俄罗斯并将其分解为所有相邻的陆地区域（POLYGONS）
Russia_exploded <- world_robinson %>%
翻译部分结束

<details>
<summary>英文:</summary>

Use `lwgeom::st_split()` as suggested in the comments, in the form of `lwgeom::st_split(POLYGON, LINESTRING)`. One thing to bear in mind is that the borders you are using are quite detailed, hence depending on the resolution of `Russia.sf`(e.g. how detailed are the borders) it may produce undesired artifacts, specially in the near the Caspian Sea (where the limits are expected to match the landmass borders).

See a full example using `giscoR` world map with a reasonable good resolution (I tested with `rworldmap` and that&#39;s when artifacts appear):

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

# Using giscoR instead
world_sf &lt;- 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)

<!-- -->

^{Created on 2023-06-28 with reprex v2.0.2}