Skip to content

proper rhumb segments #15

@mdsumner

Description

@mdsumner

Here's a decent equ-distant logic for segmentizing tiles. Underneath is a segment-decomposition of an extent, so this could be applied generally in silicate as well.

ll_extent <- function(ext, ndiscr = 24, mindist = 1e5) {
  ## technically we don't need meridionally segmentation, but this is general enough 
  ## for any set of segments assumed to be rhumb lines so it also works for graticules and
  ## will be fun for silicate
  dat2 <- matrix(c(spex::xlim(ext), 
                   spex::ylim(ext))[c(1, 2, 2, 2, 2, 1, 1, 1, 
                                      3, 3, 3, 4, 4, 4, 4, 3)], 
                 ncol = 2)
  l <- purrr::map(seq(1, nrow(dat2), by = 2), 
                  ~{
                    xy <- dat2[c(.x, .x + 1), ]
                    dst <- geosphere::distRhumb(xy[1, ], xy[2, ])
                    nn <- if (is.null(mindist)) ndiscr else round(dst/mindist)
                    nn <- max(c(nn, 3))  ## there's got to be a limit
                    cbind(approx(xy[,1], n = nn)$y, approx(xy[,2], n = nn)$y)                
                  })
   do.call(rbind, lapply(l, head, -1))
}
  

## segmentized with respect to distance for longlat segments (rhumb lines)
ext <- raster::extent(-40, 40, -60, -40)

## if we specify 
pts <- ll_extent(ext)
bind_1 <- function(x) rbind(x, x[1, , drop = FALSE])
library(sf)  
prj <- "+proj=laea +lat_0=-50 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
#prj <- "+proj=stere +lat_0=-90 +lon_0=0 +lat_ts=-70 +datum=WGS84"
sx <- st_transform(st_sfc(st_polygon(list(bind_1(pts))), crs = 4326), 
                   prj)
plot(sx, reset = FALSE)
plot(st_cast(sx, "MULTIPOINT"), add = TRUE)     


  

The more general silicate approach can be started with

sc <- silicate::SC0(spex::spex(ext))
idx <- do.call(rbind, lapply(sc$object$topology_, as.matrix))
dat2 <- as.matrix(sc$vertex[c("x_", "y_")])[t(idx), ]

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions