s2 for spherical geometry

[edzer]

Google developed a library for spherical geometry, s2, the C++ implementation of which is no longer maintained. Ege @rubak developed an R package s2 with an R interface to this library; the package contains the whole library, so has no new external dependency. Ege's UseR! talk in Brussels is found here. A nice explanation of s2, and in particular its ultra-clever indexing mechanism, is found here.

Over the weekend, I worked on a first prototype for interfacing sf with s2 (R code, tests). So far, we would get for a polar polygon wrong intersections and centroids:

> library(sf)
Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.2, lwgeom 2.3.3 r15473
> p0 = rbind(c(0,80), c(120,80), c(240,80), c(0,80))
> p1 = rbind(c(0,85), c(120,85), c(240,85), c(0,85))
> 
> pol0 = st_sfc(st_polygon(list(p0)), crs = 4326)
> pol1 = st_sfc(st_polygon(list(p1)), crs = 4326)
> st_intersection(pol0, pol1)
although coordinates are longitude/latitude, it is assumed that they are planar
Geometry set for 0 features 
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
> st_centroid(pol0)
Geometry set for 1 feature 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 120 ymin: 80 xmax: 120 ymax: 80
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POINT (120 80)
Warning message:
In st_centroid.sfc(pol0) :
  st_centroid does not give correct centroids for longitude/latitude data

With s2, we get correct values:

> sf:::s2_intersection(pol0, pol1)
Geometry set for 1 feature 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -120 ymin: 85 xmax: 120 ymax: 85
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POLYGON ((-120 85, 0 85, 120 85, -120 85))
> st_centroid(pol0)
Geometry set for 1 feature 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0 ymin: 90 xmax: 0 ymax: 90
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POINT (0 90)

The question is how to go from here. Current issues:

• s2 doesn't adopt DE-9IM or simple features; polygon intersections are only polygons, if two polygons touch this is not returned (or counted) as an intersection; there is a contains, but not the whole list of simple feature predicates
• the sf/s2 intersects/intersection code is now incredibly much slower than the (indexed) sf/GEOS implementation; not unlikely that I haven't figured out how to benefit from s2's indexing
• GEOS computations on longlat data now always warn or inform about violated assumptions; spherical geometry is the way to go to do this right
• sf already uses package geosphere for several spherical problems (but not intersects/intersections etc)
• are there alternatives to s2 we should look into? (e.g. PostGIS geography, but does that reside in liblwgeom?)
• should spherical geometry be made optional, when is which performance loss acceptable?

Reactions welcome!

[tim-salabim]

How does the perfiormance loss of s2 compare to that of geosphere?

[edzer]

geosphere does not provide intersects/intersection.

[rubak]

• I agree that it is unfortunate that s2 doesn't adopt DE-9IM or simple features, so in terms of documentation we probably have to emphasise this so users know that it is different than the rest of sf.

• I can look into using the indexing feature for testing intersections. This should hopefully speed that part up.

• I looked around for a spherical geometry library a year or two ago when I bumped into s2 for the first time. I really didn't find many alternatives.

[tim-salabim]

I'm sure you've already looked into boost?
http://www.boost.org/doc/libs/1_65_1/libs/geometry/doc/html/index.html

[edzer]

Did you find any indication that it does spherical geometry?

[rubak]

I didn't find spherical calculations in my old searches, but that might have changed or I have overlooked it. This issue mentions spherical area, so maybe there is something in there...

[rubak]

Some more details may be found here:
https://archive.fosdem.org/2017/schedule/event/geo_boost_geography/attachments/slides/1748/export/events/attachments/geo_boost_geography/slides/1748/FOSDEM17_vissarion.pdf

[tim-salabim]

I didn't do a thorough search through the documentation, but this seems to suggest that it does calculations on the sphere
http://www.boost.org/doc/libs/1_65_1/libs/geometry/doc/html/geometry/quickstart.html#geometry.quickstart.non_cartesian

also here it mentions that area can be calculated on the sphere
http://www.boost.org/doc/libs/1_65_1/libs/geometry/doc/html/geometry/reference/algorithms/area/area_1.html#geometry.reference.algorithms.area.area_1.behavior

[edzer]

You're right - thanks!

[tim-salabim]

Boost.Geometry is included in bh
https://github.com/eddelbuettel/bh/tree/master/inst/include/boost/geometry

[rubak]

I guess the big question is whether they do spherical intersections and unions of polygons. I couldn't find an example right off the bat. What do you think?
http://www.boost.org/doc/libs/1_65_1/libs/geometry/doc/html/geometry/reference/algorithms/intersection/intersection_3.html

[tim-salabim]

In the sym_difference page it says that spherical is not yet supported. I couldn't find anything for intersection
http://www.boost.org/doc/libs/1_65_1/libs/geometry/doc/html/geometry/reference/algorithms/sym_difference/sym_difference_4_with_strategy.html#geometry.reference.algorithms.sym_difference.sym_difference_4_with_strategy.behavior

[edzer]

I tried the following example:

#include <iostream>
#include <deque>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>

#include <boost/foreach.hpp>
namespace bg = boost::geometry;

int main()
{
	typedef bg::model::polygon<bg::model::point<float, 2, bg::cs::spherical_equatorial<bg::degree> > > pol;

    pol green, blue;

    boost::geometry::read_wkt(
        "POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2,3.7 1.6,3.4 2,4.1 3,5.3 2.6,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3)"
            "(4.0 2.0, 4.2 1.4, 4.8 1.9, 4.4 2.2, 4.0 2.0))", green);

    boost::geometry::read_wkt(
        "POLYGON((4.0 -0.5 , 3.5 1.0 , 2.0 1.5 , 3.5 2.0 , 4.0 3.5 , 4.5 2.0 , 6.0 1.5 , 4.5 1.0 , 4.0 -0.5))", blue);

    std::deque<pol> output;
    boost::geometry::intersection(green, blue, output);

    int i = 0;
    std::cout << "green && blue:" << std::endl;
    BOOST_FOREACH(pol const& p, output)
    {
        std::cout << i++ << ": " << boost::geometry::area(p) << std::endl;
    }
    return 0;
}

modified from their intersection example; it seemed to work OK, but I didn't compare to sf/s2. I couldn't find out well whether polygons on s2 are/can be indexed.

[edzer]

Hi @barendgehrels, long time no see! We're evaluating adopting boost::geometry (b::g) for use in R package for simple features (here) and elsewhere. Work so far focused on GEOS and google's s2. A few short questions related to b::g:

• does b::g as released (1.56.0) allow for intersections with spherical geometry?
• does b::g as released allow for indexed geometries with spherical geometries?
• does b::g as released come with WKB readers and writers?

Looks like great work!

[barendgehrels]

Hi @edzer nice to hear from you! Thanks for your interest in Boost.Geometry

You probably mean as released (1.65)?

We are currently busy with spherical geometry and geography. It is currently not me doing that, but Adam and Vissarion. Spherical segments can be intersected but the work is not all finished, so two spherical polygons cannot yet be intersected. However, that is to be expected.

WKB is (actually already a long time) provided as an extension and never released officially (contrary to WKT which was there from the beginning). I will discuss when we will release this.

About indexed geometries, @awulkiew does know all the details of this, can you comment?

So not all finished. But it would be cool to see Boost.Geometry used in R!

[awulkiew]

Hi @edzer @barendgehrels,

All relational and set operations, is_valid, is_simple, envelope, distance etc. works with spherical and geographic geometries. Aside from one specific issue that needs fixing in 1.65 (polygons containing poles) the support is quite complete. Trying it out and sharing your thoughts would definietly help us.

For geographic geometries you'd have to pass the strategies explicitly into the algorithm but they are not documented for now. By default spherical strategies are used also for geographic geometries so the results are the same in both coordinate systems.

Regarding the spatial indexing, I'm not sure what do you mean by "indexed geometries". As it is now the rtree requires the user to create envelopes, the spatial index, query the index for envelopes and then check corresponding geometries. Non-cartesian coordinate systems are supported by the rtree since it relies on relational operations of "simple" geometries like Point/Box or Box/Box. In order to use it with other geometries like Polygon you also need spherical envelope and relational operations for "more complex" geometries and as I mentioned above the support is there.

[rubak]

Definitely sounds like we should have an eye on Boost.Geometry. I have been playing a bit with using indexed queries for intersections with the s2 package. This does speed up calculations a bit for large problems, but not in any way comparable to the speed of GEOS.

I guess spherical geometry and making an index on the sphere is just more complicated than in 2D... What kind of envelope and index is used in Boost.Geometry @awulkiew ? Is everything bounding boxes and R-trees in 3-dimensional Euclidean space? Or is there a special index on the sphere?

[awulkiew]

@rubak In Boost.Geometry envelopes of 2D spherical geometries are 2D MBRs on the surface of the sphere being cartesian products in the spherical coordinates space, so they are not 3D cartesian (ECEF) boxes. They are MBRs having edges going along parallels and meridians, they are not polygons. Therefore in spherical case v.s. cartesian:

1. the envelope calculation is more complex since the max latitude of a spherical segment has to be calculated
2. spatial relations of boxes takes into account wrapping of longitudes but besides that it's close to 2D cartesian

[rubak]

@awulkiew Thanks for clarifying this. If I understand you correctly:

1. It is expected that calculations are somewhat slower in the spherical case (more complex MBRs), but fundamentally in Boost.Geometry a modified R-tree (taking wrapping into account) is used, so finding pairs of intersecting geometries should scale roughly like an R-tree search?

2. Everything is ready in Boost.Geometry to do this, so it is up to the user to test it out?

[awulkiew]

@edzer Do I understand correctly that your points are not corresponding to my points? ;)

1. Yes, the types of operations that are performed in cartesian and spherical are the same. In spherical they are only more demanding. The dimensionality of the problem and complexity of algorithms stays the same.

2. Yes, and some users are already using it.

[edzer]

Ah, yes, so it seems b::g takes care of the dateline but otherwise is identical to the case where long/lat are taken as 2D? With spherical geometry, I had in mind what S2 does: assume that edges between all vertices are great circle segments. This would mean that POLYGON((0 0,60 0,60 60,0 60,0 0)) would include POINT(30,61), which for s2 it does, but boost::geometry it doesn't? (I tried but couldn't get the simplest examples to work). Do I understand correctly that POLYGON((0 80,120 80,240 80,0 80)) with cs::spherical_equatorial degree coordinates has zero area in boost::geometry?

[rubak]

@edzer Maybe your zero area example is just related to the comment by @awulkiew about a bug for polygons containing the pole? As I understand it b::g does spherical geometry like S2 with polygon sides defined by great circle segments, but for the R-tree index it uses (lon,lat)-boxes which are not spherical polygons.

[edzer]

This is what I tried:

#include <iostream>
#include <deque>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>

#include <boost/foreach.hpp>
namespace bg = boost::geometry;

int main()
{
    typedef bg::model::polygon<bg::model::point<float, 2, bg::cs::spherical_equatorial<bg::degree> > > pol;
    typedef bg::model::point<float, 2, bg::cs::spherical_equatorial<bg::degree> > point;

    pol green, red;

    boost::geometry::read_wkt("POLYGON((0 0,60 0,60 60,0 60,0 0))", green);
    boost::geometry::read_wkt("POLYGON((30 61,31 61,31 62,30 61))", red);

    std::deque<pol> output;
    boost::geometry::intersection(green, red, output);

    std::cout << boost::geometry::wkt(output[0]) << std::endl;
    return 0;
}

but I get completely lost in the error messages from g++ 5.4.0 and boost 1.65.

[awulkiew]

In Boost.Geometry Boxes/MBRs/envelopes are not Polygons. Polygon's edges indeed are great-circle arcs however Box's edges are meridian and parallel (small-circle) arcs.

For instance:

bg::model::box<point> b;
bg::envelope(bg::model::segment<point>{{0, 60},{60, 60}}, b);
std::cout << bg::dsv(b) << std::endl;
bg::envelope(bg::model::linestring<point>{{0, 0},{0, 60},{60, 60},{60, 0},{0, 0}}, b);
std::cout << bg::dsv(b) << std::endl;

generates:

((0, 60), (60, 63.4349479675))
((0, 0), (60, 63.4349479675))

So these are Boxes with edges along parallels and meridians capable to contain the geometry. Unfortunately it seems that there is an error in spherical envelope() for areal geometries, i.e.:

bg::envelope(green, b);
std::cout << bg::dsv(b) << std::endl;
bg::envelope(red, b);
std::cout << bg::dsv(b) << std::endl;

generates:

((0, 0), (60, 60))
((30, 61), (31, 62))

which is wrong. The result should be the same as above.

Regarding the compilation error, you're passing std::deque<> which was not adapted to any Boost.Geometry concept (like MultiPolygon) into bg::wkt(). So either use bg::model::multi_polygon<pol> or pass output[0] into the bg::wkt().

Regarding the question about the area of POLYGON((0 80,120 80,240 80,0 80)) this polygons has coordinate in invalid range (240 should be -120, though for area this is not a problem but in general this is not correct) and is counter-clockwise (but your type is clockwise). For modified: POLYGON((0 80,-120 80,120 80,0 80)) I get: 0.0399250108917.

[barendgehrels]

@awulkiew thanks for comments and explanations!

[edzer]

Thanks for the explanation. I updated the example above, and made both rings counter clockwise. Still, the intersection I get is an empty polygon, which surprised me.

[awulkiew]

Yes, the result should not be empty. For me the following program (tested Boost 1.65.0 and 1.65.1 with msvc-14 and mingw-w64 g++4.8):

typedef bg::model::point<float, 2, bg::cs::spherical_equatorial<bg::degree> > point;
typedef bg::model::polygon<point, false> pol;

pol green, red;
boost::geometry::read_wkt("POLYGON((0 0,60 0,60 60,0 60,0 0))", green);
boost::geometry::read_wkt("POLYGON((30 61,31 61,31 62,30 61))", red);

std::deque<pol> output;
boost::geometry::intersection(green, red, output);

std::cout << boost::geometry::wkt(output[0]) << std::endl;

generates:

POLYGON((30 61,31 61,31 62,30 61))

Is it different for you?

[rubak]

It seems like Boost.Geometry is probably the way to go for spherical calculations in sf. It seems like a very good fit in the sense that it knows about WKT/WKB and adopts DE-9IM (as far as I can tell). My C++ skills are not good enough to try to understand a new C++ library, so I will leave the implementation to someone else and stick with s2 for my own purposes until Boost.Geometry is adequately interfaced from R. I think this is really exciting and hope you can find someone with a solid C++ background to interface Boost.Geometry from R.

[edzer]

@rubak I agree, although we haven't explored CGAL, which I believe now also comes as header templates.

@awulkiew it turns out I compiled against boost 1.58 (ubuntu xenial); I had in mind 1.65 because that is the BH R package with boost headers, which R packages link against (I mostly use c++ for R packages, which pick up 1.65, barely "rare" as above). My apologies!

@barendgehrels let me know when you want me to look into WKB I/O, and where I can find it.

[rubak]

@edzer I didn't know that CGAL also is a header library now. That is very interesting. I have looked at it several times, but it has always been too big a system dependency for me. Maybe that drawback is over now. However, after a bit of searching it appears that the default way of installing CGAL even as header only uses CMAKE and linking to gmp, mpfr, and zlib is recommended, so there might be some work in making this work with the standard R package toolchain. Anyway, it's another interesting option for someone to look deeper into...

[awulkiew]

@edzer I'm glad everything works as expected.

Regarding WKB, all extensions has to be grabbed from Boost.Geometry develop branch. In this case from here:
https://github.com/boostorg/geometry/tree/develop/include/boost/geometry/extensions/gis/io/wkb

[edzer]

I moved the s2 stuff to a branch for now. I will look into the boost stuff deeper when I understand template programming better -- in sf we handle geometries at the geometry level, not at the level of their types. I can see the examples for this in the boost docs, I just can't read them, yet.

[barendgehrels]

OK @edzer , thanks for trying it out. Actually the intention of Boost.Geometry is that users do not need any knowlegde of the underlying templates or design. Apart from defining the type (maybe we could also avoid that), no templates are necessary for basic operations.

But I understand, as soon as you get errors, it is often hard to figure out the real cause between all template noise. Especially if you use the brand new features from the library.

Anyway, if you need more help, just ping.

Regards, Barend

[SymbolixAU]

If anyone's interested I too wanted to learn a bit more about Boost Geometry, so decided to start implementing it in R in my repo.