Points/polygons/polylines should be reverse-projected

[jcheng5]

From the tutorial materials here:

http://www.uclm.es/profesorado/vgomez/useR2015/

BMcD <- readOGR("datasets/BMcD.shp", "BMcD")
leaflet(BMcD) %>% addTiles() %>% addMarkers()

This ought to work but doesn't, due to BMcD using a non-WGS84 projection. The CRS info is included on the BMcD object, though, so we should do the conversion automatically. Here is a manual way:

coords <- rgdal::project(BMcD@coords, BMcD@proj4string@projargs, inv = TRUE)

[edzer]

Since rgdal::project does not take care of datum transformation, and leaflet (I assume) will need WGS84 coordinates, and a more robust approach is to use spTransform:

require(sp)
BMcD <- spTransform(BMcD, CRS("+ellps=WGS84 +proj=longlat +datum=WGS84 +no_defs"))

but still, this puts several of the locations in the river (where they are not!) The problem (not with leaflet) seems to be that the useR2015 material from Virgilio contains a shapefile with a wrong CRS; the following would resolve this:

library(leaflet)
library(sp) # for meuse, but also for spTransform
demo(meuse, ask=FALSE, echo = FALSE) # loads meuse

BMcD <- readOGR("datasets/BMcD.shp", "BMcD")
proj4string(BMcD) = proj4string(meuse) # overwrite CRS
BMcD <- spTransform(BMcD, 
  CRS("+ellps=WGS84 +proj=longlat +datum=WGS84 +no_defs"))

leaflet(BMcD) %>% addTiles() %>% addMarkers()

which puts the points on the river banks.

[rsbivand]

The BMcD shapefile has the following WKT:

"+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 
+x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"

which was from a version of the EPSG database current about ten years ago, when the files were created. The "meuse" dataset in sp has:

"+init=epsg:28992 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel 
+towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs"

which crucially has a +towgs84= tag, and the most recent EPSG version included in PROJ.4 4.9.1 has:

+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 
+x_0=155000 +y_0=463000 +ellps=bessel 
+towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,
-1.87740163998045,4.0725 +units=m +no_defs

that is 28992 with added precision in the +towgs84 terms.

In relation to the original question, please note that WGS84 is the World Geodetic System of 1984, defining an ellipsoid and a general datum. It has been revised on a number of occasions since publication. It is however NOT a projection, but part of a definition of a geographical coordinate reference system (usually in decimal degrees). Projected coordinates, such as the Meuse Bank data, need the definition of the projection in addition to that of the underlying geographical coordinates.

With reference to the comment on the original posting, one should NEVER project and/or transform coordinates automatically unless one is completely sure about the provenance and authority of the definitions of the source and target spatial reference systems, or willing to accept the ensuing uncontrolled inaccuracy.

[jcheng5]

one should NEVER project and/or transform coordinates automatically

@rsbivand Given that Leaflet needs the end-result to be in EPSG:3857 (Web Mercator) in order to visualize it, what would you recommend instead of converting automatically? Throwing an error requesting that the user please project to EPSG:3857 first? Do the automatic conversion but emit a warning?

[becarioprecario]

Hi,

Many thanks for your comments. I have simplified Edzer's code to use epsg:28992:

library(leaflet)
library(sp)
library(rgdal)

BMcD <- readOGR("datasets", "BMcD")
proj4string(BMcD) <- CRS("+init=epsg:28992")
BMcD <- spTransform(BMcD, CRS("+ellps=WGS84 +proj=longlat +datum=WGS84 +no_defs"))

leaflet(BMcD) %>% addTiles() %>% addMarkers()

However, I have not been able to produce a shapefile with the correct CRS using writeOGR, as it seems to drop "+init=epsg:28992".

@jcheng5, you may simply show a message stating that data have been re-projected. I believe this is what plotKML() does, for example. As Roger pointed out, this may not be the best solution but at least you will tell the user that data have been re-projected.

Best,

Virgilio

[chris-holcomb]

I think something might be getting lost here, as Leaflet ingests Vector layers such as Polygon, and Geojson as lat lng's in WGS84 (EPSG:4326), not Web Mercator (EPSG:3857). Leaflet actually transforms these from EPSG:4326 to EPSG:3857 for you.

I.e., the proj4 CRS for EPSG:3857 is actually +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs (http://epsg.io/3857). This is sometimes also wrongly referred to 900913. WGS84 is +proj=longlat +datum=WGS84 +no_defs (http://epsg.io/4326).

In Leaflet itself you could skip this processing step by simply showing the map in EPSG:4326. The other options are EPSG:3395 and something called "Simple", besides that you can use to Proj4Leaflet for other projections. A very nice trick would be to change the projection to the best state plane upon selection or zooming to a specific state or county. More info: http://leafletjs.com/reference.html#icrs.

However, Leaflet will not let you mix input EPSGs in your basemaps or WMS Overlays and you should use EPSG:3857 for those unless you over-ride the output projection system.

TL;DR: vector data could be re-projected to EPSG:4326 with a warning. ogrInfo pulls the original proj4 from a shapefile ala ogrInfo(dsn=system.file("vectors", package = "rgdal")[1], layer="cities")[["p4s"]]

[jcheng5]

Yeah, I understand--the question was whether the automatic wgs84 to epsg3857 projection should be called out to the user via a warning or message. (I would not have thought it necessary, but Roger and Virgilio have surely forgotten more about spatial data analysis than I will ever know...)

[chris-holcomb]

Automatic EPSG:X to EPSG:4326 should show a warning (it's not terribly reliable), but there's never a need to project into EPSG:3857 (on the R side).

[jcheng5]

I know my code doesn't project into EPSG:3857, but sending it to Leaflet with the default settings effectively is performing the projection.

But I like your rule... @rsbivand @becarioprecario would you agree with that? Converting to EPSG:4326 (if necessary) will show a warning, but then the data will be sent to Leaflet to be displayed in EPSG:3857 without any further message/warning?

[chris-holcomb]

Gotcha, yea I think that falls into Leaflet's purview and I haven't seen any issues going from 4326 to 3857 (it's getting it out of some obscure projection that causes issues). If you are going this far, something else to consider in the transformation is optional simplification (default FALSE) via gSimplify, as this can greatly reduce file size. I also restrict the number of digits in my lat lng's but I'm not sure off-hand how to do this in R (I'm pulling from PostGIS). BTW, absolutely tremendous work on this package; I'm finalizing our transition this week after spending forever on data prep.

[edzer]

It makes sense to convert non-WGS84 longlat Spatial* objects to epsg 4326 before sending them to leaflet, and warn in that case.

You can find out whether an object x deriving from Spatial is epsg:4326 by the test

!is.projected(x) && grep("+datum=WGS84", x@proj4string@projargs)

which should be TRUE; other possible outcomes are FALSE (conversion needed) and NA (CRS unknown).

[rsbivand]

Sorry to appear infrequently, I'm largely offline to mid-August. The real question remains the provenance and authority of EPSG:X. Very often the data do not declare the complete set of information required, especially the +towgs84= coefficients or the +datum=. This means that although we know the target CRS, we don't know the input CRS sufficiently accurately to proceed. In some jurisdictions the input CRS may in fact be restricted - the Indian military keep some details confidential. It is important to make the user aware of the issue and responsible for resolving any uncertainties, as this cannot be done reliably in an automatic way. If the background contextual map has a known CRS, and this is the target, using a guessed or incomplete input CRS will yield a mismatch, either giving wrong positional information, extracting wrong altitude or other values, putting points in the sea/lake, and so on. Some input data even mistakes Eastings and Northings. The key issue is then not something that can be coded (ever). It is to alert the user to their need to understand that incomplete input CRS only permit plotting without contextual information, that is that their "map" should never be combined with any other map until the CRS is made sufficiently complete.

By the way, "+datum=WGS84" is not always present in a correct proj4 declaration, as it may be rendered as "+towgs84=0,0,0,0,0,0,0" and "+ellps=WGS84", so this cannot be used for checking. Note that users will often try to correct CRS wrongly, or in the US setting guess the wrong State Plane CRS; in Europe the ED50 datum has many +towgs84= versions (as with the Meuse Bank data set) - they change in precision and authority over time.

[rsbivand]

An addendum: I'm minded to provide a better look-up function than just using grep() on the output of rgdal::make_EPSG() - which uses the PROJ.4 epsg.csv file. Ideas?

A further possibility is to automate access to the appropriate Grids & Datums column at:

http://www.asprs.org/Grids-Datums.html

where these are relevant for the target CRS:

http://www.asprs.org/a/resources/grids/02-2014-World.pdf http://www.asprs.org/a/publications/pers/2015journals/PERS_May_2015/HTML/index.html#356

(the latter no longer easy to read). The idea would be to solicit a country code, or maybe catch a click on a world map, and return a link to the G&D column(s), Again, ideas?

Unfortunately, I think that spatialreference.org is not actively maintained.

[bhaskarvk]

proj4leaflet is now merged in master, it would be unadvisable to auto transform CRS to EPSG4326, It is best for users to do it on their own before passing data to leaflet. This is an added extra step but better done intentionally than doing it internally w/ a warning.

More documentation at https://github.com/bhaskarvk/leaflet/issues/16#issuecomment-249369147.

I'm closing this issue but feel free to open again or add more thoughts.