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 r-spatial/sf#1649 (comment)

[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.

[Nowosad]

@mtennekes we spoke with @edzer about this issue and we think it would be great to have this for the tmap4 release.

Also -- similar work done for interactive maps -- https://kalberer.org/pirmin/blog/equal-area-projection/

[mtennekes]

@Nowosad: do you mean the idea at the end of this discussion, so an automatic crs recommendation based on the area (in case the crs of spatial object is in WGS84)? Yes, definitely! What would be a good place for it? sf, tmap, tmaptools, ...?

Also it would still be nice to draw orthogonal globes like in the initial post and the one from @edzer https://r-spatial.org/book/08-Plotting.html#sec-transform. For that we'll need to check if that code works for other orthogonal projections as well. I mean with the aim to "cut the earth from space" (#929).

[edzer]

There's also s2::s2_plot():

ne = rnaturalearth::ne_countries(50)
s2::s2_plot(st_make_valid(ne))

from @paleolimbot

[edzer]

[Nowosad]

@Nowosad: do you mean the idea at the end of this discussion, so an automatic crs recommendation based on the area (in case the crs of spatial object is in WGS84)? Yes, definitely! What would be a good place for it? sf, tmap, tmaptools, ...?

Not even a recommendation, but a default. Quoting @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." So, when the user wants to visualize data in, e.g., EPSG:4326, on a global scale then the output map will use Eckhart IV by default, but if the data in EPSG:4326 is just for some region -- then an orthographic projection would be used...

[tim-salabim]

FWIW, I will investigate integrating https://github.com/maplibre/maplibre-gl-leaflet into leafem which should support https://kalberer.org/pirmin/blog/equal-area-projection/
However, I am not sure at this stage how many baselayers are available for this projection at low zoom levels. But, with vector tile baselayers, we may be able to reproject on the fly using some jsvascript implementation of PROJ?
@edzer may I ask, why "Eckhart IV for global maps" is a better choice than equalearth?

[mtennekes]

Yes, agree with having that as a default.

Once we have this function, say optimal_crs(x), how can we use this in tmap @Nowosad ? An option is to set tmap_options(crs = "auto").

@tim-salabim My thoughts exactly about the availability of basemaps. I'm open to any good pseudo cylindrical crs. See #929 (comment) for the most common ones. Perhaps eqearth?

[tim-salabim]

Sorry, @mtennekes I meant equalearth, not naturalearth... Fixed in previous comment, so 👍🏽

[Nowosad]

Once we have this function, say optimal_crs(x), how can we use this in tmap @Nowosad ? An option is to set tmap_options(crs = "auto").

👍🏻

[mtennekes]

I've started a new branch, autoCRS.

See https://github.com/r-tmap/tmap/blob/autoCRS/R/auto_crs.R for the general auto crs function. Happy to migrate to another package (sf, tmaptools) if that makes more sense.

A simple heuristic:

IF Is the area of the bounding box at least 1/4 of the earth surface? AND Is the mean latitude between -50 and 50?
THEN pseudo cylindrical crs
ELSE orthogonal

of course open to improvements....

Apart from the function itself, important is that is works well in tmap, so please test it. The intended behaviour:

• option crs set to "auto" for plot mode, meaning that any shape for which sf::st_is_longlat() will be set to this auto crs.
• the user can set crs = "auto" in tm_shape to force it to be drawn in this auto projection, even if the shape itself is projected.
• (to be implemented) the user should also be able to set crs = "original" meaning that the crs of the shape is used, and not an automatic crs, even if the shape is unprojected.

Some code to test:

# projected shape, so using the original crs (epsg 28992)
tm_shape(NLD_prov) +
	tm_polygons() +
	tm_graticules()

# force auto crs
tm_shape(NLD_prov, crs = "auto") +
	tm_polygons() +
	tm_graticules()

tm_shape(World) +
	tm_polygons() +
	tm_graticules()

# only working with s2=TRUE (with warnings) due to graticules
tm_shape(World[World$continent == "Africa", ]) +
	tm_polygons() +
	tm_graticules()

tm_shape(World[World$continent == "South America", ]) +
	tm_polygons() +
	tm_graticules()

tm_shape(World[World$name == "Norway", ]) +
	tm_polygons() +
	tm_graticules()

tm_shape(World[World$name == "China", ]) +
	tm_polygons() +
	tm_graticules()


# below only working well with s2=FALSE

northern_hami = sf::st_bbox(c(xmin = -180, xmax = 180, ymin = 0, ymax = 90), crs = 4326) |> sf::st_as_sfc()
World_north = sf::st_intersection(World, northern_hami)

southern_hami = sf::st_bbox(c(xmin = -180, xmax = 180, ymin = -90, ymax = 0), crs = 4326) |> sf::st_as_sfc()
World_south = sf::st_intersection(World, southern_hami)

tm_shape(World_north) + 
	tm_polygons() +
	tm_graticules()

tm_shape(World_south) + 
	tm_polygons() +
	tm_graticules()

tm_shape(northern_hami) + 
	tm_polygons() +
	tm_graticules()

tm_shape(southern_hami) + 
	tm_polygons() +
	tm_graticules()

[edzer]

IF Is the area of the bounding box at least 1/4 of the earth surface? AND Is the mean latitude between -50 and 50?
THEN pseudo cylindrical crs
ELSE orthogonal

this will probably fail with areas that are small but strongly N-S elongated.

[edzer]

I looked at your auto_crs proposal; some more observations/suggestions:

My idea of using an auto crs was strictly for unprojected data: if data are projected, you'd use that projection. The point is, unprojected data need to be projected to plot them, projected data not.

Do we want to cater for areas crossing the antemeridian? In that case current bbox is not useful (will result in mean_lon=0). s2 has a function that returns a more useful bbox:

library(s2)
s2::s2_bounds_rect(s2_data_countries("Antarctica"))
#   lng_lo lat_lo lng_hi    lat_hi
# 1   -180    -90    180 -63.27066
s2::s2_bounds_rect(s2_data_countries("Netherlands"))
#     lng_lo   lat_lo   lng_hi  lat_hi
# 1 3.314971 50.80372 7.092053 53.5104
s2::s2_bounds_rect(s2_data_countries("Fiji"))
#    lng_lo    lat_lo    lng_hi    lat_hi
# 1 177.285 -18.28799 -179.7933 -16.02088

where lng_lo > lng_hi indicates straddling the antemeridian.

When, for an orthographic projection the center is close to one of the poles, you'd probably want to switch to one of the polar projections.

[mtennekes]

My idea of using an auto crs was strictly for unprojected data: if data are projected, you'd use that projection. The point is, unprojected data need to be projected to plot them, projected data not.

We're on the same page: in tmap, auto_crs is only used to project unprojected data.

Thx for the s2_bounds_rect tip! Very useful.

Additional ideas:

• I also thought about calculating the lat and long range distances in km, and deriving the aspect ratio of the spatial object/bbox as a proxy to choose the CRS. Of course this is only useful when the poles are not inside the bbox.
• Yes, indeed better to snap to polar projections if the center is close to one of the poles. We could also snap to other standard/official CRSs, e.g. if bbox ~ Europe, use LAEA (EPSG 3035). However, not sure though where to stop if we go down that road...

Do you have suggestions @tim-salabim @mdsumner? By the way, Michael: what CRSs are commonly used in Australia for world maps? The standard ones, or also ones with Australia in the center?

[tim-salabim]

My 2 cents regarding web mapping:

• if you look at the leaflet example of combinded web mercator and equalearth projections, things don't really behave nicely when the jump between the projections occurs (zoom into greenland for example) - https://equal.bbox.earth/leaflet-eq2merc/
• we would need some heuristic for reprojecting all overlay data along with the basemaps if we go down this route... Though I think the code is already there.

In my mind, the best solution I've seen so far is what google maps does - gradually projecting to a rotating globe, the further you zoom out. Of course this will prohibit views where you want to see all data overlaid over the entire globe... All in all, I am not sure what the best solution is for web mapping...

[mdsumner]

We use the usual ones, and there's a mix of local UTM, regional and local LCC same as other places have.

Interesting to advise Orthographic for small areas, that's fine and is what Manifold.net used to have everywhere in the docs, I think it's better to use either LAEA or AEQD or ORTH and make it explicit why you would choose that - you can zoom our pretty far with reasonable result everywhere and not notice the difference, but with LAEA and AEQD you at least have a valid reason for choosing (exact area - within topology limits of your geometric representation of course - or exact distance from anywhere to the centre point). For elongated meridional regions I use OMERC, but that's harder to automate because the params are different.

I'd definitely insert LCC for mid-sized regions, but you need extra heuristics to choose the latitude range secants. People always say STERE for polar, but you can orient them anywhere and LAEA is a good choice too if you want actual area rather than (local) shape preserved.

Please don't confine this to tmap - put the logic somewhere reusable? You could do this with points and base plot and PROJ, show why choices were made, scope for improvement and expansion, and then port to general format usage with wk::wk_transform, and fold it into here for the narrower cases.

You could also use the warper to demonstrate, say from a standard global topo like GEBCO (500m), and even patch in COP30m when you're really local, it's a lot more robust than vector handling in R will for arbitrary conversions.

[mdsumner]

also paging @h-a-graham who's been looking at family-local reprojection helpers

in SOmap we have an _auto function that takes a proj family (the short string laea, aeqd, laea, aea, stere, etc) and fits it to the input data, I always wanted to break that out into a minimal general thing, and let the user play around with the family properties

https://github.com/AustralianAntarcticDivision/SOmap/blob/main/R/SOmap_auto.R

[mtennekes]

Thx for your input @mdsumner, very useful!

Agree that this should not belong totmap (or even tmaptools) but at a lower level. But where?

• sf
• PROJ
• other package, namely: ...
• new (lightweight) package

[h-a-graham]

Thanks for the tag @mdsumner. This sounds like a great idea and has a lot of use more broadly, from my perspective - here's a gist I put together which is relevant to this discussion: https://gist.github.com/h-a-graham/3e20fb61dd1f1287ae25a748e4581e1e

For context, my intention with this was to be able to simply select an appropriate CRS to build satellite image composites anywhere in the world using a sensible projection.

FWIW, I think this could fit really nicely into a tiny new package which could have two main exported functions, something akin to what is in the above gist, providing a bit more control and then an auto_crs wrapper which contains the scaling logic for making the appropriate choice based on a general rule set as discussed elsewhere in this discussion. These functions could just return the crs allowing it to be agnostic of the different rspatial classes; other packages could then transform as they wish with as few or as many additional proj inputs as desired. 😄

[mtennekes]

Back to the (original) topic: orthogonal projections now implemented. Based on Edzer's code, but had to a bit more to make it work for other centers. Code still ugly and possibly unstable (?): https://github.com/r-tmap/tmap/blob/master/R/misc_crs.R

library(tmap)

plot_ortho = function(lat, lon) {
    tm_shape(World, crs= paste0("+proj=ortho +lat_0=", lat, " +lon_0=", lon), bbox = "FULL") +
        tm_polygons("HPI", fill.scale = tm_scale_intervals(values = "cols4all.pu_gn_div", midpoint = 27.5)) +
        tm_style("natural")
}

plot_ortho(lat = 30, lon = 0)   

plot_ortho(lat = 30, lon = 120) 

plot_ortho(lat = -10, lon = -60)    

plot_ortho(lat = 90, lon = 0)   

plot_ortho(lat = -90, lon = 0)

^{Created on 2024-11-04 with reprex v2.1.0}

Now this projection is also used for step 3 (transformation), which is undesirable. Better to use this crs just for the plotting phase (step 4). I'll open a new issue for this.

[mtennekes]

Just discovered that this only works with the option earth.boundary = TRUE (which is enabled in the style "natural"). Should work without too.

tm_shape(World, crs = "+proj=ortho +lat_0=0 +lon_0=0") +
    tm_polygons()