Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active March 14, 2019 05:13
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mdsumner/2be160dcd66f35893e4e7c926f264c22 to your computer and use it in GitHub Desktop.
Save mdsumner/2be160dcd66f35893e4e7c926f264c22 to your computer and use it in GitHub Desktop.
st_recenter <- st_recentre <- function(x, clon = NULL, ..., tryfix = TRUE) {
if (is.null(clon)) return(x)
if (!st_is_longlat(x)) stop("recentring not appropriate for non longlat data")
## try to fix problematic geometry
if (tryfix) {
if (all(grepl("POLYGON", st_geometry_type(x)))) x <- suppressWarnings(st_buffer(sf::st_as_sf(x), 0))
x <- st_wrap_dateline(x, options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
}
wbox <- st_as_sfc(st_bbox(c(xmin = -180, ymin = -90, xmax = (clon)%%360 - 180, ymax = 90), crs = st_crs(x)))
west <- suppressWarnings(st_intersection(x, wbox))
west <- st_set_geometry(west, st_geometry(west) + c(360, 0))
west <- st_set_crs(west, st_crs(x))
east <- suppressWarnings(st_crop(x, st_bbox(c(xmin = (clon)%%360 - 180, xmax = 180, ymin = -90, ymax = 90))))
xx <- rbind(
west, east
)
## ensure geometries are of consistent type
xx <- sf::st_cast(xx)
xx
}
@mdsumner
Copy link
Author

Examples

  st_recenter <- st_recentre <- function(x, clon = NULL, ..., tryfix = TRUE) {
  if (is.null(clon)) return(x)
  if (!st_is_longlat(x)) stop("recentring not appropriate for non longlat data")
  ## try to fix problematic geometry
  if (tryfix) {
   if (all(grepl("POLYGON", st_geometry_type(x)))) x <- suppressWarnings(st_buffer(sf::st_as_sf(x), 0))
   x <- st_wrap_dateline(x, options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
  }
  wbox <- st_as_sfc(st_bbox(c(xmin = -180, ymin = -90, xmax = (clon)%%360 - 180, ymax = 90), crs = st_crs(x)))
  west <- suppressWarnings(st_intersection(x, wbox))
  west <- st_set_geometry(west, st_geometry(west) + c(360, 0))
  west <- st_set_crs(west, st_crs(x))
  east <- suppressWarnings(st_crop(x, st_bbox(c(xmin = (clon)%%360 - 180, xmax = 180, ymin = -90, ymax = 90))))
  xx <- rbind(
    west, east
 ) 
  ## ensure geometries are of consistent type
  sf::st_cast(xx)
  xx
}


## other scales don't work, but this dataset is LARGE
# x <- rnaturalearth::ne_countries(returnclass = "sf", scale = 10)

## LINES
x <- rnaturalearth::ne_coastline(returnclass = "sf")
library(sf)
#> Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0
plot(st_recentre(x, -120))
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning in st_is_longlat(x): bounding box has potentially an invalid value
#> range for longlat data

#> Warning in st_is_longlat(x): bounding box has potentially an invalid value
#> range for longlat data

## POINTS
plot(st_recentre(st_cast(x, "MULTIPOINT"), -120))
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning in st_is_longlat(x): bounding box has potentially an invalid value
#> range for longlat data

#> Warning in st_is_longlat(x): bounding box has potentially an invalid value
#> range for longlat data

## POLYGONS
data("wrld_simpl", package = "maptools")
plot(st_recentre(st_as_sf(wrld_simpl)[7], 147))
#> dist is assumed to be in decimal degrees (arc_degrees).
#> Warning in CPL_wrap_dateline(st_geometry(x), options, quiet): GDAL Error
#> 1: IllegalArgumentException: Points of LinearRing do not form a closed
#> linestring
#> Warning in CPL_wrap_dateline(st_geometry(x), options, quiet): GDAL Error
#> 1: IllegalArgumentException: Points of LinearRing do not form a closed
#> linestring
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar

library(ggplot2)
ggplot(st_recentre(st_as_sf(wrld_simpl), 147)) + geom_sf(aes(fill = POP2005, colour = NAME)) + xlim(75, 250) + guides(colour = FALSE)
#> dist is assumed to be in decimal degrees (arc_degrees).
#> Warning in CPL_wrap_dateline(st_geometry(x), options, quiet): GDAL Error
#> 1: IllegalArgumentException: Points of LinearRing do not form a closed
#> linestring

#> Warning in CPL_wrap_dateline(st_geometry(x), options, quiet): GDAL Error
#> 1: IllegalArgumentException: Points of LinearRing do not form a closed
#> linestring
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar

Created on 2019-03-14 by the reprex package (v0.2.1)

@mdsumner
Copy link
Author

A few improvements

  st_recenter <- st_recentre <- function(x, clon = NULL, ..., tryfix = TRUE) {
  if (is.null(clon)) return(x)
  if (!st_is_longlat(x)) stop("recentring not appropriate for non longlat data")
  ## save the crs while we do our munging
  crs <- st_crs(x)
  x <- st_set_crs(x, NA)
  
  
  ## try to fix problematic geometry
  if (tryfix) {
   if (all(grepl("POLYGON", st_geometry_type(x)))) x <- suppressWarnings(st_buffer(sf::st_as_sf(x), 0))
   x <- st_wrap_dateline(x, options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
  }
  wbox <- st_bbox(c(xmin = -180, ymin = -90, xmax = (clon)%%360 - 180, ymax = 90))
  west <- suppressWarnings(st_crop(x, wbox))
  west <- st_set_geometry(west, st_geometry(west) + c(360, 0))
  east <- suppressWarnings(st_crop(x, st_bbox(c(xmin = (clon)%%360 - 180, xmax = 180, ymin = -90, ymax = 90))))
  xx <- rbind(
    west, east
 ) 
  ## ensure geometries are of consistent type
  xx <- sf::st_cast(xx)
  bb <- st_bbox(xx)
  ## hmmm
  # if (bb["xmax"] > 180 && !grepl("\\+over", crs) && !grepl("init", crs)) {
  #   crs <- sprintf("%s +over", crs)
  # }
  st_set_crs(xx, crs)
}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment