Tmap: Orthogonal projections

[mtennekes]

Playing around with orthogonal projections...

library(tmap)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

# function to check if polygons contain 4 points
sf_is_valid_4poly = function(x) {
    g <- sf::st_geometry(x)
    vapply(g, function(gi) {
        nrow(st_coordinates(gi)) == 5L
    }, FUN.VALUE = logical(1))
}

# function to create world surface
world_surface = function(datum = 4326, step = 2, nx = 360/step, ny = 180/step, projection = NULL, union = TRUE) {
    s = stars::st_as_stars(sf::st_bbox(c(xmin=-180,ymin=-90,xmax=180,ymax=90), crs = datum), nx = nx, ny = ny)
    s2 = sf::st_as_sf(s)
    s4 = if (!is.null(projection)) {
        s3 = sf::st_transform(s2, crs=projection)
        s3[sf_is_valid_4poly(s3), ]
    } else{
        s2
    }
    if (union) {
        sf::st_union(sf::st_buffer(s4, dist = 1e-03))
    } else {
        s4
    }
}

# transform World and create surface
data(World)
ortho_crs = st_crs("+proj=ortho +lon_0=0 +lat_0=35")
World_ortho = st_transform(World, crs = ortho_crs)
World_ortho <- st_make_valid(World_ortho)
surface = world_surface(projection = ortho_crs) #takes 10-15 seconds

tm_shape(surface) +
    tm_polygons("lightskyblue1") +
    tm_shape(World_ortho) + 
    tm_polygons("darkolivegreen3") +
    tm_graticules(x = seq(-180, 150, by = 30), y = seq(-90, 90, by = 30), labels.show = FALSE) +
    tm_layout(frame = FALSE)
#> Warning: The shape World_ortho contains empty units.

^{Created on 2020-06-07 by the reprex package (v0.3.0)}

A couple of things for which I need your ideas and help:

• Some orthogonal crs's (with other lon_0 and lat_0 values) gave errors. This one was error-free, but not without problems; as you can see, China and the US disappeared. (@edzer: any tips?)
• world_surface is slow, because it transforms a medium sized stars object to an sf. Any ideas how to improve it? So what I want to achieve, is a large polygon that covers the whole earth. I wrote a function called tmaptools::bb_earth a few years ago, but this doesn't work with orthogonal projections.
• Finally, it would be nice to have an interactive globe as well (in view mode). Perhaps threejs::globejs could be used. Other suggestions (or pull requests:-)) are welcome!

[edzer]

Yes, see https://github.com/r-spatial/sf/blob/master/demo/twitter.R for how messy this can get. I wouldn't go this way until @paleolimbot and I get libs2 done; see https://github.com/paleolimbot/s2plot - I expect this to be operational in a few months. If you want to get here faster, consider helping us with libs2 ;-)