leaflet does not handle projected coordinates and confuses naive users

[rsbivand]

A user asked me why sp and rgdal are broken because leaflet does not display sp objects read in using rgdal. It is because leaflet() simply does not explain in the help pages (there is no vignette) how to display a completely valid object (cause: not in web mercator - a completely broken idiocy - or longlat). After digging on the off-CRAN hard to navigate github site, I found addPolygons(), and transformed the data. Provide proper guidance, please.

library(rgdal)
dsn <- system.file("vectors", package = "rgdal")[1]
ogrInfo(dsn=dsn, layer="scot_BNG")
scot_BNG <- readOGR(dsn=dsn, layer="scot_BNG")
library(leaflet)
addPolygons(leaflet(scot_BNG))
addPolygons(leaflet(spTransform(scot_BNG, CRS("+init=epsg:4326"))))
library(mapview)
mapview(scot_BNG)

[mdsumner]

Where is this documented in sp, for the behaviour of adding to a plot paying no heed to the coordinate system previously used? That would be a good place to find inspiration. I'm happy to PR this doc item.

[jcheng5]

Thanks for the feedback, @rsbivand. (cc @schloerke)

This has come up multiple times:

#139
#349
r-spatial/mapview#72

Each time there has been strong opinions on several sides but it looks like the latest discussion on mapview ended with everyone happy. @tim-salabim it looks like mapview now calls spTransform for all sp objects with a CRS (that isn't an exact match for 4326) and st_transform for all sf objects. Any reason for Leaflet not to do the same at this point?

[tim-salabim]

The current mapview implentation calls st_transform(x, 4326) on all x with non-NA crs. sp objects are converted to sf before that step.

[rsbivand]

Thanks for considering this. I think there are three choices: an informative warning for non-conformant input; use proj4s inside leaflet if leaflet supports it; use sf or rgdal on the R side with a large loaded image downside. The warning conditions would be needed anyway for b and c, so I'd choose a. However, I've always felt that users benefit from pushback when they input data that are non-conformant in these kinds of settings, and many stats books on spatial data have nice sections on coordinate reference systems (Waller & Gotway 2004, Banerjee et al. 2004). Once the user grasps the issue, the problem resolves itself from there.

[rsbivand]

@tim-salabim, in mapview/R/projection.R, you had checks for string identity to check and warn, but don't use try() or similar. Not only does proj4 4.* have issues with some transformation to llcrs, but they will get more pronounced with proj4 5.*, which implements a geodetic pipeline without the WGS84 datum hub. I think that just going to 4326 is OK, but for 5.0.0 a workaround is in place in rgdal for web mercator (proj4 5.0.0 correctly places 3D data up to 10km below or above the WGS84 surface!).

I see that while leaflet imports raster, raster only suggests rgdal and sf, so keeping the load footprint sensible. I think it is best for leaflet to either throw an error or warn on NA crs and non-conformant crs on input (and then only raster, sp and sf input, telling the user to transform first), rather than importing the GDAL i/o parts of sf and rgdal that are not needed.

[tim-salabim]

@rsbivand these checks were only in place for a short time in the dev version. The discussion at r-spatial/mapview#72 suggested that the easiest is to simply and blindly transform to 4326 whatever comes in as there is no overhead for data already in 4326. If there are issues re proj4 5.* could you provide an example where transformation doesn't work as expected in the mapview repo?
As I have said before, I agree that leaflet does not need to provide on-the-fly reprojection (though I have no objections if it would) but an informative warning would be great.

[rsbivand]

All we have on proj.4 5.0.0 is that web mercator was broken, and we have reports on grass-dev and gdal-dev of differences. I cannot handle all the issues that will arise, the update in proj4 was overdue, and WGS84 is increasingly not fit for purpose because the world has changed in shape a good deal since then. I'm watching GDAL and GRASS, see r-spatial/sf#545.

[rsbivand]

For this issue, it would be ideal if leaflet threw an informative error on data that are non-conformant. I do not like the mapview approach as a wrong crs will throw errors that the user is not expecting. I'd certainly expect to see sf::st_transform() and rgdal::spTransform() (not sp::) inside a try() or tryCatch() to give an error message that will alert the user to the problem (no or erroneous crs that cannot be transformed). There was a lot of this with geom_sf, where inverse projections that did not exist were attempted to draw graticules (IMO, it shouldn't draw graticules, waste of ink).

[jcheng5]

Thanks all for your thoughts. I talked to @schloerke and we've decided that since there are some big decisions to make here, we're going to release the current master branch to CRAN as-is, but definitely tackle the issue in a minor release following soon after. We're so many months late in updating to Leaflet.js 1.x that we don't want to delay it any more, nor do we want to rush to make a decision without fully understanding all the tradeoffs.

[rsbivand]

Good, that makes a lot of sense. @tim-salabim - could we try to find the trip-point criteria to establish when action should be taken? @edzer - do we know when a declared crs is "valid", and so can be transformed to a epsg:4326 or Web Mercator hub even though latest proj has abandoned fixed transformation hubs?

[tim-salabim]

@rsbivand sure, though I need a bit of guidance here. I am happy to wrap things in try or tryCatch and provide a more informative error/warning if we can identify clearly when things go wrong.
Is the issue of a "wrong" web mercator projection both ways? In mapview we transform from web mercator to 4326 if data is in web mercator. This is then fed through the leaflet package to leaflet.js which does the transformation to web mercator (at least for vector data). @jcheng5 @schloerke does leaflet.js use proj to transform to web mercator?

[rsbivand]

It was wrong in proj 5.0.0, and is trapped in rgdal (so guarded against). 5.0.1 was released 36 hours ago and the issue should be resolved there: OSGeo/PROJ#834. I'll report back when checked.

[jcheng5]

@tim-salabim Leaflet.js does its own projection /~https://github.com/Leaflet/Leaflet/tree/master/src/geo/projection

[rsbivand]

The proj 5.0.1 bug-fix release resolves the web-mercator issue. I see that Leaflet.js does not seem to use proj4 in any form unless an addin is present, so should be isolated from ongoing work on proj itself.

[Christoph999]

@rsbivand Not sure whether that fits here, because it is a way lower level, but I give it a try: Does that have the same root cause? I am really unable to understand what happens, although I spent quite some time with search in the www.

[rsbivand]

@Christoph999 no, not really. leaflet is for doing interactive displays with low-level access. You will find the units and spatial representation issues easier in sf I think, using mapview or tmap for interactive mapping.

> set.seed(1)
...
> library(sf)
Linking to GEOS 3.7.2, GDAL 3.0.0, PROJ 6.1.0
> library(mapview)
> obj <- st_as_sf(data.frame(lng=lng, lat=lat, id=1:n), coords=c("lng", "lat"))
> st_crs(obj) <- 4326
> mapview(obj)

If you then need to transform, use st_transform().

[Christoph999]

@rsbivand Thanks a lot for your input!