Sf: Question about processing message on st_intersects

[lbusett]

A short question:

checking for intersections between sf object in geographical coordinates always yield the message: 

"although coordinates are longitude/latitude, it is assumed that they are planar"

Can I ask why is it so ? Could this lead to "incorrect" behaviour or can it be safely suppressed/ignored ?

[rsbivand]

Depending where you are on the globe, and the distance between points forming a line or polygon, an apparently straight line in geographical coordinates will be more or less curved when projected to the plane. So if you are lucky, there may be no impact (small study area close to the equator?) or you may be unlucky. The warning is relevant.

[edzer]

Good question. Here is an example:

> arctic = st_polygon(list(rbind(c(0,80), c(120,80), c(240,85), c(0,80))))
> pole = st_point(c(0,90))
> st_intersects(pole, arctic)
[[1]]
integer(0)

Arguably, this is right, because sf has no clue what coordinates refer to - they could be local coordinates on a 2D brain scan slice.

In this example:

> st_intersects(st_sfc(pole, crs = 4326), st_sfc(arctic, crs = 4326))
although coordinates are longitude/latitude, it is assumed that they are planar
[[1]]
integer(0)

sf is informed that these are long/lat Earth coordinates, and the arctic polygon is supposed to include the pole, but st_intersects says it doesn't, which is wrong. It's wrong because it assumes coordinates are planar, and you are warned about that assumption.

(reality is somewhat more complicated: every polygon on the sphere divides the sphere into two parts, depending on the ring direction; this would have to be taken into account)

I hope to extend sf with routines that do proper intersections on the sphere; R package s2 from Ege @rubak seems to be a good candidate for this.

[lbusett]

Thanks for the answers. Would it make sense to promote the Message to a Warning (maybe with some additional explanation specifying which/where problems may occurr) ?

[edzer]

In principle yes, but for many use cases (small area not near the poles) it is of little importance; everyone with such cases would have to take care of all the warnings for all geom ops, making scrips unreadable.

[lbusett]

Ok, good. I'd suggest however to add some more specification to the message, so that users may be more "confident" on wheter they can ignore and suppress it. What about something like:

"Although coordinates are longitude/latitude, it is assumed that they are planar. This may lead to errors/inaccuracies when considering large areas or working on areas near the poles"

[edzer]

This still leaves many questions open; I think the current message exactly tells what is going on.

[Schuch666]

Hey @edzer, I have a question:

I am using a continent-size shape to distribute gases emissions on a global grid, this assumption make much difference? there is a way to estimate this relative error for this case and for a city scale?

PS: no areas near the poles

[edzer]

I don't think that your question addresses this issue, I also don't see how it relates to sf.

[jzadra]

Seems that this is still around in 2021 - is there a recommended transform prior to st_intersection that will negate this problem?

[edzer]

Maybe see https://www.r-spatial.org/r/2020/06/17/s2.html and https://r-spatial.github.io/sf/articles/sf7.html

[jzadra]

Will do, thanks!