st_transform does not transform map projection correcttly in some situations

[PanfengZhang]

Describe the bug st_transform does not transform map projection correcttly in some situations.

To Reproduce

library(sf)
library(ggplot2)
library(rnaturalearth)
wld <- ne_coastline(returnclass = "sf")
wld2 <- st_transform(wld, crs = "+proj=natearth")
ggplot() + 
  geom_sf(data = wld2)

When I use the function of st_transform in sf package to transform map projection from "wgs1984" to "Natural Earth"(https://proj.org/operations/projections/natearth.html), the horizontal line appear on the the high latitudes of the northern hemisphere. The figure shown as follow.

library(sf)
library(ggplot2)
library(rnaturalearth)
wld <- ne_countries(returnclass = "sf")
wld2 <- st_transform(wld, crs = "+proj=natearth +lon_0=110")
ggplot() + 
  geom_sf(data = wld2) +
  theme_bw()

When I transform the Longitude of projection center, I'm getting more inexplicable errors. The figure shown as follow.

Why is that?

[rsbivand]

The lines are not broken at the edge (dateline) longitudes, so on projection wrap around, the problem is not with st_transform() but the input geometries:

wld3 <- st_transform(st_wrap_dateline(wld), crs = "+proj=natearth")
ggplot() +
  geom_sf(data = wld3)

The documentation of st_wrap_dateline() mistakenly shows "DATELINEOFFSET=180" as a way to solve your second problem, but it isn't. This option only shows the tolerance in +/- degrees around -180 and 180 degrees, default +/-10, see https://gdal.org/programs/ogr2ogr.html#cmdoption-ogr2ogr-wrapdateline et seq. That is, my comment and example in https://github.com/r-spatial/sf/issues/541 seems not to hold any longer.

[PanfengZhang]

Thank you

[PanfengZhang]

library(rnaturalearth)
library(sf)
library(ggplot2)
wld <- ne_coastline(returnclass = "sf") |> st_wrap_dateline() |> 
  st_transform(crs = "+proj=moll +lon_0=120")
ggplot(data = wld) + geom_sf()

If the world map changes the projection center, there will still be an error, is there a good solution to solve this issue now?

[rsbivand]

It appears that a custom intersection between the input and a split bounding box is a possible resolution:

library(sf)
library(ggplot2)
library(rnaturalearth)
wld <- ne_coastline(returnclass = "sf")
st_crs(wld) <- "OGC:CRS84"
wld2 <- st_transform(wld, crs = "+proj=natearth +lon_0=110")
ggplot() + 
  geom_sf(data = wld2) +
  theme_bw() # problem
bb0 <- st_bbox(wld)
bb1w <- bb0
bb1e <- bb0
bb1w[3] <- -70 - 0.0001 # antimeridian of target minus fuzz
bb1e[1] <- -70 + 0.0001 # plus fuzz
bb1 <- c(st_as_sfc(bb1w), st_as_sfc(bb1e))
sf_use_s2(FALSE)
wld_70W <- st_intersection(wld, bb1) # avoid s2 because we need a planar bounding box
wld2 <- st_transform(wld_70W, crs = "+proj=natearth +lon_0=110")
ggplot() + 
  geom_sf(data = wld2) +
  theme_bw()

[PanfengZhang]

Thank you. It solved my problem.

[zimmerberg-ms]

Segmentation of the bounding boxes produces smooth borders (relevant for e.g. Mollweide and Robinson projections)

bb1 <- bb1 %>% st_set_crs(NA) %>% st_segmentize(1) %>% st_set_crs(4326)

This would also be a great addition for st_break_antimeridian()