Dissolve internal boundaries of multipolygon that crosses 180th meridian

Question

I am trying to make a map using the world dataset from spData using a polar projection. Everything looks great, except that there are internal boundaries within Russia where it crosses the 180th meridian, and it appears that st_union() won't dissolve polygons across that boundary, at least as I am using it. How can I dissolve the boundaries so that they don't show up on my map? I would prefer answers that use sf.

library(sf)
library(spData)
library(dplyr)
library(ggplot2)

# polygon to crop the countries to only the northernmost ones
crop_poly <- tibble(geometry = st_sfc(st_point(c(0, 90)),
                                      crs = 'EPSG:4326')) %>%
  st_sf() %>%
  st_transform(crs = 'EPSG:3413') %>%
  st_buffer(dist = 3909167) %>%
  smoothr::densify(n = 3) %>%
  st_transform(crs = 'EPSG:4326')

# crop the world dataset
world_north <- world %>%
  st_intersection(crop_poly) %>%
  st_transform(crs = 'EPSG:3413') %>%
  select(name_long, continent, region_un, subregion) %>%
  st_union(by_feature = TRUE) # this is the suggested method for dissolving internal boundaries

# make the map
ggplot(world_north) +
  geom_sf(data = crop_poly,
          color = 'transparent',
          fill = 'gray90',
          alpha = 0.5) +
  geom_sf(color = 'black', fill = 'gray80') +
  geom_sf(data = crop_poly, 
          color = 'gray70',
          fill = 'transparent',
          inherit.aes = FALSE) +
# add an arrow to highlight location of undesired internal boundaries
  geom_segment(x = 0, xend = -1300000,
               y = 0, yend = 1300000,
               arrow = arrow(),
               color = 'red') +
  theme_bw()

Other things I have tried

I have also tried playing around with st_wrap_dateline(), although I'm not really sure if this is an application in which it could be useful, but I couldn't figure out a way that it would solve the issue. I attempted to use suggestions from this post to solve my problem, but I saw no difference (the internal boundaries did not get dissolved) when I tried various combinations of EPSG:32635, st_union(), and st_simplify(). The various permutations of code I tried are not included here to keep this post readable.

Well, first of all, that is **NOT** the international date line, as you would have known if you'd bothered to look. It does appear to be the 180th longitude, so try looking at what longitude lines do or don't do in `sf` . — Carl Witthoft, Aug 25 '23 at 18:41
@CarlWitthoft, thank you for pointing out that the terminology I used wasn't the most accurate. I have edited the question. I would like you to know, however, that I found your comment to be very disrespectful to me and the time and effort that I put into writing this question. — Heidi Rodenhizer, Aug 25 '23 at 19:59
You should maybe try the GIS stack exchange site for R-spatial questions: https://gis.stackexchange.com/ — Spacedman, Aug 26 '23 at 08:16
@Spacedman, yeah, that's a good point. I have written a few questions over there, but I have less reputation, which can be rather annoying. — Heidi Rodenhizer, Aug 28 '23 at 21:55

score 5 · Accepted Answer · answered Aug 26 '23 at 07:58

Russia is the fourth feature in your data, lets pull that out for convenience:

russia = world_north$geom[4]

Its a MULTIPOLYGON:

> russia
Geometry set for 1 feature 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -2235484 ymin: -1622414 xmax: 3909167 ymax: 3909167
Projected CRS: WGS 84 / NSIDC Sea Ice Polar Stereographic North
MULTIPOLYGON (((-1448716 1400198, -1507001 1384...
>

So lets split it to look at the individual polygons:

u = st_cast(russia,"POLYGON")

and (by various inspections, eg eventually plot(u[c(5,12)])) we can get the distance between the two parts across the 180 lat line:

> st_distance(u[c(5,12)])
Units: [m]
          [,1]      [,2]
[1,] 0.0000000 0.4020763
[2,] 0.4020763 0.0000000

Telling us there's a 40.2cm gap between them, so they won't be dissolved. That could be a very important piece of 40.2cm wide territory. Or a numerical precision problem.

If you don't mind risking an international incident, you can add 25cm to the Russian boundary all-round and then the line disappears as sf reconstructs the MULTIPOLYGON geometry to be valid - no union needed (the solution with less geopolitical ramifications is to edit the polygon at the vertex level, just changing some points along the join - easy to do interactively with a GIS package):

plot(russia)
plot(st_buffer(russia,0.25))

This has also taken out the 180 border in the small island, I suspect these parts were also centimetres apart.

You can fix this directly into the data frame:

world_north$geom[4] = st_buffer(world_north$geom[4],.25)

and we're done:

that small island is the [Wrangel Island](https://en.wikipedia.org/wiki/Wrangel_Island), the final place on Earth to support woolly mammoths as an isolated population — Jindra Lacko, Aug 28 '23 at 16:20
For further context, the slight offset between polygons which cross the 180th meridian was a conscious choice to avoid geometry issues with s2. See [this issue](https://github.com/Nowosad/spData/issues/67#issuecomment-1699502457) for the details. — Heidi Rodenhizer, Aug 30 '23 at 17:05

Dissolve internal boundaries of multipolygon that crosses 180th meridian

Other things I have tried

1 Answers1