Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Created November 26, 2021 23:38
Show Gist options
  • Save mdsumner/f181204c4c3b7f09aa4a2111e07f07a8 to your computer and use it in GitHub Desktop.
Save mdsumner/f181204c4c3b7f09aa4a2111e07f07a8 to your computer and use it in GitHub Desktop.
  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)

@mdsumner
Copy link
Author

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

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