the cylinder -to sphere thing

[mdsumner]

Just a code dump for now tocylr does direct translation from planar (igh in this case) to cylinder

  list(cbind(x, y)[c(seq_along(x), 1), ])
}
cut <- sf::st_sfc(sf::st_multipolygon(list(
  to_polygon(rep(c(80.01, 79.99), each = 91),   c(-90:0, 0:-90)), # third cut bottom
  to_polygon(rep(c(-19.99, -20.01), each = 91),   c(-90:0, 0:-90)), # second cut bottom)  
  to_polygon(rep(c(-99.99, -100.01), each = 91),   c(-90:0, 0:-90)), # first cut bottom)
  to_polygon(rep(c(-40.01, -39.99), each = 91), c(90:0, 0:90) )
)
), crs = 4326)

data("wrld_simpl", package = "maptools")
## we have to buffer, else topology problems (other sources are not perfect either)
world_sf <- sf::st_buffer(sf::st_as_sf(wrld_simpl), 0)
library(sf)
x <- x <- st_difference(world_sf$geometry,cut)
#plot(st_transform(x, "+proj=igh"))     

gr <- graticule::graticule(seq(-180, 180, by = 10), seq(-90, 90, by = 10), nverts = 120, tiles = TRUE)
grc <- st_difference(st_as_sf(gr), cut)
#plot(st_transform(grc, "+proj=igh"))

## avoid sf preventing overplot
plot(st_geometry(st_transform(grc, "+proj=igh")), col = rgb(0.95, 0.95, 0.95, 0.5), border = "transparent")
plot(st_cast(st_transform(grc, "+proj=igh"), "LINESTRING"), add = TRUE, col = rgb(0.65, 0.65, 0.65))
plot(st_transform(world_sf, "+proj=igh"), add = TRUE, col = "transparent", border = "black")

## now mesh and plot in rgl
library(anglr)
im <- ceramic::cc_location(cbind(0, 0), buffer = 6378137 * pi)
mesh0 <- anglr::as.mesh3d(DEL0(st_transform(x, "+proj=igh"), max_area = 1e10), image_texture = im)
#mesh_plot(mesh)

mesh <- mesh0
tocyl <- function(lonlatheight, rad = 6378137) {
  theta = lonlatheight[,1] * pi/180
  v = lonlatheight[,2] * pi/180
  x = rad * cos(theta)
  y = rad * sin(theta)
  cbind(x, y, v)
}
tocylr <- function(lonlatheight, rad = 6378137) {
  theta = lonlatheight[,1] 
  v = lonlatheight[,2] 
  x = rad * cos(theta)
  y = rad * sin(theta)
  cbind(x, y, v * rad/pi)
}

scl <- function(x) scales::rescale(x, to = c(-pi, pi))
#xyz <- reproj::reproj(t(mesh$vb[1:3, ]), source = "+proj=igh", target = 4326)
xyz <- cbind(scl(mesh$vb[1, ]), scl(mesh$vb[2, ]), scl(mesh$vb[3, ]))
xyz[!is.finite(xyz)] <- NA
mesh$vb[1:3, ] <- t(tocylr(xyz))

gmesh <- SC0(st_cast(st_transform(grc, "+proj=igh"), "LINESTRING"))
xyz <- cbind(scl(gmesh$vertex$x_), scl(gmesh$vertex$y_), 0)
xyz[!is.finite(xyz)] <- NA
gmesh$vertex$z_ <- 0
gmesh$vertex[c("x_", "y_", "z_")] <- tocylr(xyz)
glb <- anglr::globe(SC0(world_sf))
plot3d(mesh)
plot3d(gmesh, add = TRUE)
plot3d(glb, add = TRUE, color = "lightgrey")
aspect3d("iso")
bg3d("black")

play3d(spin3d(c(0, 1, 1)))

https://twitter.com/mdsumner/status/1295246118087880704?s=20

https://twitter.com/mdsumner/status/1295154142969647104?s=20