sf
sf copied to clipboard
st_transform does not transform map projection correcttly in some situations
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?
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.
Thank you
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?
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()
Thank you. It solved my problem.
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()