proper rhumb segments

[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), ]

[Maschette]

How does this work if you want multiple polygons such as #14 ?

[mdsumner]

We'd have to loop it over all tiles

g <-graticule::graticule(lons = c(95, 135, 175),lats = c(-42, -39), tiles = TRUE)

bind_1 <- function(x) rbind(x, x[1, , drop = FALSE])

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))
}

x <- sf::st_sfc(purrr::map(split(g, 1:nrow(g)), ~sf::st_polygon(list(bind_1(ll_extent(.x))))), crs = 4326)

plot(st_transform(x, projection(Bathy)))

It could be neater, but I can't handle it rn

[Maschette]

Cheers, is three any reason tiles=TRUE seems to return multiple polygons now? Is this a new thing, I dont remember coming across is when we were developing SOmap

library(graticule)
#> Loading required package: sp
g <-graticule::graticule(lons = c(-180.0, -124.0 ), lats = c(-42, -40), tiles = TRUE)
plot(g)

^{Created on 2019-10-04 by the reprex package (v0.3.0)}

[mdsumner]

something dumb on my part, forget graticule and create a raster with the exact extent and grain you want:

r <- raster::raster(extent(95, 175, -42, -39), ncols = 2, nrows= 1, crs = "+init=espg:4326")
g <- spex::qm_rasterToPolygons_sp(r)

then as before with g