st_as_sf faithfully uses "+proj=longlat +ellps=WGS84" instead of conversion to "+proj=longlat +datum=WGS84"

[charliejhadley]

Currently I use geosphere to sample great circles as I can't figure out how to do it with sf directly, which looks as follows:

library("geosphere")
library("sf")
library("leaflet")
sp_thing <- gcIntermediate(
  p1 = c(13.96013, 50.97178),
  p2 = c(-74.00594, 40.71278),
  n = 50,
  addStartEnd = TRUE,
  sp = TRUE
)
sf_thing <- st_as_sf(sp_thing)
leaflet() %>%
  addPolylines(data = sf_thing)
## Warning message:
## sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +no_defs).
## Need '+proj=longlat +datum=WGS84' 

Should st_as_sf faithfully take

sp_thing@proj4string
## CRS arguments: +proj=longlat +ellps=WGS84 

Or should it be clever and coerce it to "+proj=longlat +datum=WGS84" silently?

Please note - I recognise this is a little bit of a Rube Goldberg machine and that this issue might be closed as such, if so thanks for reading anyhow

[edzer]

A problem with leaflet, as

> sf_thing
Simple feature collection with 1 feature and 0 fields
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: -74.00594 ymin: 40.71278 xmax: 13.96013 ymax: 55.94954
epsg (SRID):    4326
proj4string:    +proj=longlat +ellps=WGS84 +no_defs
                        geometry
1 LINESTRING(13.96013 50.9717...

has a perfectly valid WGS84 proj4string. @tim-salabim is this known in the r-leaflet community?

Alternatively, you could use

mapview::mapView(sf_thing)

which does not complain.

All three proj4strings in your post point to the same coordinate reference system; doing literal string matches to one of them doesn't help; the order of the arguments are for instance allowed to change.

[charliejhadley]

Thanks for comprehensive reply, appreciated. I'll post an issue at /~https://github.com/rstudio/leaflet about this and close this one.

I like mapview! But I've got a bunch of leaflet templates I can't move away from right now.

Cheers,

Martin

[edzer]

FYI (see also here):

ibrary("sf")
sf_thing <- st_linestring(rbind(
  c(13.96013, 50.97178),
  c(-74.00594, 40.71278))) %>%
  st_sfc(crs = 4326) %>%
  st_segmentize(units::set_units(50, km))

also creates your sf_thing; behind the scenes it calls geosphere functions.

[tim-salabim]

Well, this (or somewhat related things) have surfaced over the past in both mapview and leaflet. I am unsure whether the leaflet folks are aware of the problem of literal vs. semantic. In any case, I just added a comment regarding this in the above mentioned leaflet issue.
It is a tricky thing to handle in leaflet as it is not resptricted to data of certain spatial classes (as opposed to mapview) so cannot by default take care of projections. But emitting warnings where no warning is needed is also not ideal.

[edzer]

A sufficient test would be to find a match both to "+proj=longlat" and to "WGS84" somewhere in the string.

[tim-salabim]

That's what I suggested:
r-spatial/mapview@e97c7af#diff-5798b2525d12114a21aa684a5bb163f8R120

[edzer]

Looks more complicated that what I proposed: why for instance insist on the +nodefs?