Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Meridian-wrap artefacts in shapefile reprojected with sf::st_transform() #1495

Closed
geotheory opened this issue Sep 24, 2020 · 3 comments
Closed

Comments

@geotheory
Copy link

I'm getting 180° meridian-related artefacts when projecting a shapefile to Robinson - see horizontal red lines in the map below. Any suggestions why this is occuring or how to resolve?

require(sf)
require(ggplot2)

m = read_sf('https://gist.githubusercontent.com/geotheory/b75a60e1659daa484e0f513a6350be79/raw/3838e475af1cdfa73bd94f68fd0807ed3e9c956f/robinson-artefacts.geojson')
m = st_transform(m, crs = '+proj=robin +ellps=WGS84 +datum=WGS84 +no_defs')

ggplot(m, aes(fill = value)) + 
  geom_sf(col = 'red', size=.5) +
  theme_light() + theme(panel.border = element_blank())

image

@edzer
Copy link
Member

edzer commented Sep 24, 2020

Right after reading the geojson, I see

> st_bbox(m)
      xmin       ymin       xmax       ymax 
-179.99958  -55.98125  180.00042   83.62792 

so there seems to be a problem with your input data.

@geotheory
Copy link
Author

geotheory commented Sep 24, 2020

Good spot @edzer. Do you know the best practice way to rectify this? I find that the (hacky)

for(i in 1:length(m$geometry)) m$geometry[[i]] = m$geometry[[i]] + st_point(c(-0.00042,0))` 

does in practice move the polygons, but st_bbox(m) does not update.

@mdsumner
Copy link
Member

Try crop the input to -180/180 after also zero-buffering (that decomposes the shapes to segments and forces a degeneracy fix for some problems, and appears to work with your data).

Hacky? Sure, but balance that with what you need to do ;)

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(ggplot2)

m = read_sf('https://gist.githubusercontent.com/geotheory/b75a60e1659daa484e0f513a6350be79/raw/3838e475af1cdfa73bd94f68fd0807ed3e9c956f/robinson-artefacts.geojson')
m = st_buffer(m, 0)  ## fixes some issues
#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
#> endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).
m = st_transform(st_crop(m, st_bbox(c(xmin = -180, xmax = 180, ymin = -90, ymax = 90))),  crs = '+proj=robin +ellps=WGS84 +datum=WGS84 +no_defs')
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

ggplot(m, aes(fill = value)) + 
  geom_sf(col = 'red', size=.5) +
  theme_light() + theme(panel.border = element_blank())

Created on 2020-10-21 by the reprex package (v0.3.0)

@edzer edzer closed this as completed Mar 3, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants