Improve crs stuff

[mtennekes]

Triggered by #457

Historically in tmap, projected spatial objects were preferred, not just for plotting, but also for calculating distances and areas. (With the important message to use a proper CRS for the map area). Nowadays, calculation of distances and areas is easy with s2. Not sure about the accuracy, but at least good enough and at least better than an unsuitable CRS, and therefore safer to use.

We need to decide where to do the CRS transformation in the process: before step 3 (transformation, e.g. cartograms, centroids, etc), or before step 4 (plotting). In the current v4 implementation, it takes place before step 3, but if the crs is specified as option, it is subsequently applied in step 4 (needed for the view mode). If the keep this approach, then we need to think about how to specify the crs and bbox for step 4 only.

# this is good: crs used for geometric distortion in step 3, and also for step 4 (plotting)
tm_shape(World, crs = "+proj=robin") +
tm_cartogram(size = "pop_est")

# this is not so good, because there may be step-3 layer transformations,
# for which an orthogonal projection is not good for global use cases
tm_shape(World, crs = "+proj=ortho +lat_0=0 +lon_0=0") +
tm_polygons() +
tm_style("natural")

# better to use this crs only in step 4:
tm_shape(World) +
tm_polygons() +
tm_style("natural") +
tm_options( crs = "+proj=ortho +lat_0=0 +lon_0=0")

But this is not very user-friendly I guess. Should be add another element for this purpose @Nowosad @nickbearman and others?

E.g.

tm_crs("+proj=ortho +lat_0=0 +lon_0=0")

Which sets the crs for the whole plot, so users don't have to figure other which shape is the main shape.
If we do it this way, a follow-up: what to do with the bounding box? tm_bbox ? Or a combined element tm_plot(crs = ..., bbox = ...)

Also I/we need to improve the spatial manipulation. Probably we can make use of s2 much more. Also tmaptools::bbis still used a lot in tmap, but it is outdated so needs an update.

[mtennekes]

Current status of the implementation (work-in-progress)

Two internal variables:

• crs_step3: The crs for step 3 (transformation: think cartograms)
• crs_step4: The crs for step 4 (plotting)

Leaflet requires vectorized objects to be in 4326 and raster images in 3857. I'll make a new tmap option for that, e.g. crs_data which is not supposed to be changed by the user.

However, leaflet can still plot in projected crs's:

epsg9311 <- leafletCRS(
	crsClass = "L.Proj.CRS",
	code = "EPSG:9311",
	proj4def = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs",
	resolutions = 2 ^ (16:7)
)

leaflet(options = leaflet::leafletOptions(crs = epsg9311)) |> 
	addPolygons(data = World)

so crs_step4 is still relevant.

My proposal is that if there are no basemaps, and the shape is projected (or a crs is specified), this crs is used as crs_step4.

The tmap option crs is used to specify crs_step4, so the crs for the plot. Probably I'll create a shortcut tm_crs, because I'd like to keep tm_options for advanced users only, and tm_layout (which offers a subset of options) for styling only.

Ideally, I'd like to have a 'smart' auto crs in case a non-projected shape is used.
@h-a-graham and @mdsumner : any follow-ups on #457 (comment)

If that gist is ok, but not ready for a CRAN package yet, I'm also happy to put it into tmap for the time being. It's probably still beter in the vast majority of use cases than lat/longs.

[mtennekes]

Also: I'm looking for R-leaflet code that uses a projected-crs tile servers. The ones I had seem to be taken offline.

[tim-salabim]

Ideally, I'd like to have a 'smart' auto crs in case a non-projected shape is used.

mapview uses L.CRS.Simple for this case

This seems to work nicely for whatever projection the coordinates have. Two caveats of this approach are that

• we cannot ensure that multiple layers are plotted in a sensible manner (but that's expected and a user issue)
• we cannot draw rasterImages this way (maybe that could work with leafem::addGeotiff() which I have not tested so far.

world9311 = sf::st_transform(World, "EPSG:9311")
mapview::mapview(world9311, native.crs = TRUE)

sf::st_crs(world9311) = NA
sf::st_crs(World) = NA

mapview::mapview(world9311, native.crs = TRUE)
mapview::mapview(World, native.crs = TRUE)

storms = sf::st_as_sf(atlStorms2005)
sf::st_crs(storms) = NA

# storms located in the northern US somewhere
mapview::mapview(world9311, native.crs = TRUE) +
  mapview::mapview(storms, native.crs = TRUE)

# storms located correctly with regard to world as both have same crs coincidentially
mapview::mapview(World, native.crs = TRUE) +
  mapview::mapview(storms, native.crs = TRUE)

[h-a-graham]

Sadly I've not had any time to package up those ideas, I think if that solution works for tmap as it stands then bringing that in (or something similar) makes a lot of sense, it would be straight forward to swap out for a package when it appears... I like the idea of tm_crs, I wonder if that could have some kind of optional proj_type argument with a limited set of options that match up with the generic projection function so that a user could optionally specify the projection in a more general sense, enabling the testing of a range of appropriate projections without having to think about proj4 strings etc?

[mdsumner]

I do think it's worth having some logic to pick a crs family for a given set of longlat data. Essentially anything is fine at Tasmania size (laea, aeqd, lcc, ortho etc) and at Australia size probably lcc, then you've got special cases like polar, hemisphere, global, and very long and skinny spans along a great circle (tmerc, or omerc), but can we sensibly categorise those from raw data? Maybe we should also have bbox lookups for "nations" and "oceans" or known grids line epsg 3577 and the conus equivalent. I'll try a few things, mostly I'd like ways to promote attention to these ideas, I usually find it easy to advise in real cases. Having some analytical grids to see errors could be handy too (we could use mdsumner/tissot) I think a raster field as error indicator indicator is better than those little indicatrixes

[mdsumner]

and because, I suppose, the first level (local) is to have approximately equal area and equal distance and conformality across the entire map, which is trivial as long as you are centred on a small area - so anything will do, but next level out you really need a user to care about area (laea) or distance (aeqd) or shape (stereographic), or some compromise of those (conic at first, or any number of other specific customizations), and then there's "must be global"

I suppose there's value in having a threshold "this is so local, it doesn't matter what you choose", or this is so nearly global you actually have to have some destructive error unavoidably in some places - and in the middle you have an opportunity for folks to actually learn about the projections a bit. Having a think about "orthographic", I think I'd strike that off the list because of the clipping problem. At least with laea, aeqd, and stereographic they only have topology or numeric problems at the outer margins, whereas orthographic needs a pretty nuanced clipping and adaptive segmentation to really be done well.

[mtennekes]

Thx for the insights @mdsumner

tm_crs function updated. It now includes an automatic crs recommendation, with the message that it is work-in-progress and that it may eventually migrate to a new stand-alone package.

By default, unprojected shapes will still be plotted unprojected (plate carree), but for global enough shapes, a message with by printed:

tm_shape(World) + tm_polygons("HPI") 

[tip] Consider a suitable map projection, e.g. by adding `+ tm_crs("auto")`.

I opted for this behaviour instead of plotting it by default with this automatic crs, to make users fully aware of it: now they 'have to' act to get a proper crs.

The addition of tm_crs("auto"), so

tm_shape(World) + tm_polygons("HPI") + tm_crs("auto")

will use an automatic crs. Which one?

1. for world maps, crs_global is used, which is a tmap option, now set to +proj = "eqearth".
2. otherwise, there is an argument of tm_crs called "property", which can be set on "global" (to force using crs_global, "area" (equal area), "distance" (equidistant), and "shape" (conformal). The last three use the families "laea", "aeqd", and "stere" respectively.

I've included your gist @h-a-graham (many thanks) https://github.com/r-tmap/tmap/blob/master/R/tm_crs.R#L96 with some minor changes: stere family added, stars accepted, instead of centroid of the first polygon, the centroid of the bounding box is used for cent_coor.

Again, work in progress. New ideas welcome.

[mdsumner]

Wow nice, will check this out 🙏

[h-a-graham]

Looks great, glad the gist was of some use 🙂, looking forward to taking it for a spin!