graticule_masking.md

  library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1
edge0 <- function(x, y, ndiscr = 18) {
  dx <- if(length(x) > 1) diff(x) else diff(y)
  sf::st_set_crs(sf::st_segmentize(sf::st_sfc(sf::st_linestring(cbind(x, y))), dx/ndiscr), 
               "OGC:CRS84")
}
north <- function(x = c(-180, 180), y = 90, ndiscr = 18) {
 edge0(x, y, ndiscr)  
}

south <- function(x = c(-180, 180), y = -90, ndiscr = 18) {
  edge0(x, y, ndiscr)
}

west <- function(x = -180, y = c(-90, 90), ndiscr = 18) {
  edge0(x, y, ndiscr)
}
east <- function(x = 180, y = c(-90, 90), ndiscr = 18) {
  edge0(x, y, ndiscr)
}

proj <- "+proj=murd3 +lat_1=30 +lat_2=50"
plot(x <- st_transform(st_line_merge(st_cast(c(south(ndiscr = 180), west(), north(), east()), "MULTILINESTRING")), proj))

## to polygonize
st_cast(st_polygonize(st_union(x)))
#> Geometry set for 1 feature 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -21994670 ymin: -10018750 xmax: 21994670 ymax: 21721620
#> CRS:           +proj=murd3 +lat_1=30 +lat_2=50
#> POLYGON ((19718792 21721618, 19933541 21274455,...

^{Created on 2021-11-27 by the reprex package (v2.0.1)}

update, now only using sf for the final format object and reproject/plotting

the west, north, east, south thing is no ordered to avoid having to do a line merge (but we can do that in various ways)

  edge0 <- function(x, y, nverts = 6) {
  segmentize(cbind(x, y), nverts)
}
segmentize <- function(x, nverts = 6) {
  xx <- seq(x[1L, 1L], x[2L, 1L], length.out = nverts)
  cbind(xx, approxfun(x)(xx))
}
north <- function(x = c(-180, 180), y = 90, nverts = 18) {
 edge0(x, y, nverts)
}

south <- function(x = c(180, -180), y = -90, nverts = 18) {
  edge0(x, y, nverts)
}

west <- function(x = -180, y = c(-90, 90), nverts = 18) {
  edge0(y, x, nverts)[,2:1]
}
east <- function(x = 180, y = c(90, -90), nverts = 18) {
  edge0(y, x, nverts)[,2:1]
}

linearize <- function(x) {
  dupes <- duplicated(x)
  x[!dupes, ]
}

polygonate <- function(x) {
  x <- linearize(x)
 rbind(x, x[1L, , drop = FALSE])
}
to_sf <- function(x) {
  sf::st_sfc(sf::st_polygon(list(polygonate(x))), crs = "OGC:CRS84")
}
proj <- "+proj=murd3 +lat_1=30 +lat_2=50"
x <- rbind(west(), north(), east(), south(nverts = 90))
plot(sf::st_transform(to_sf(x), proj))

^{Created on 2021-11-27 by the reprex package (v2.0.1)}

see this stuff proceed at https://github.com/mdsumner/llbox