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

[mtennekes]

This is one for v4. See also https://github.com/r-spatial/sf/issues/1649#issuecomment-860819840

[tim-salabim]

well, cesium is still my preferred virtual globe implementation. But this would require major efforts to get it into a similar state of usability as e.g. leaflet or mapdeck... I've said it before, we'd need something like a GSOC or RConsortium project to implement it. I'd be happy to serve as a "mentor", but I won't have resources to actually implement it.

[mtennekes]

Agree. I still have negative resources myself.

I saw this advertisement, already from 5 years ago: https://github.com/rstats-gsoc/gsoc2016/wiki/cesium:-R-interface-to-cesium.js-virtual-globe Still seems like a fun project for young R-enthusiasts like @luukvdmeer.

[tim-salabim]

Oha, I totally forgot about that one :-)

[Nowosad]

Interesting update: https://github.com/goergen95/cesium

[mtennekes]

Yes, saw it. Would be nice candidate for a new mode :-)

[Nowosad]

There is a discussion about CRSs at the Spatial Data Science across Languages (SDSL) event (https://r-spatial.org/sdsl/). One suggestion there (from @edzer) is to try a different default when WGS84 is used (e.g., orthographic projection).

[edzer]

I was thinking orthographic for smaller areas, like having less than 10% of the globe, and Eckhart IV for global maps, maybe also for anything in between. And do the same thing for sf and tmap as default. What do you think?

[mtennekes]

Love this idea @edzer

Aren't there rules of thumb of which CRS to recommend for which scale?

My initial thoughts on this: orthographic for smaller areas definitely good, and pseudo-cylindrical like Eckert IV for global maps also. For intermediate zoom levels it also depends on the center location. The vast majority of these global map projections are centred around lat=0 & long=0. If it is too far off, say North America, there are better alternatives.

Would indeed be nice to have a common CRS-recommendation function somewhere in r-spatial.