Wrong union result in spherical coordinate system

[remster]

Consider the test case:

namespace cartesian {
  using point = boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>;
  using linestring = boost::geometry::model::linestring<point>;
  using polygon = boost::geometry::model::polygon<point>;
  using multipolygon = boost::geometry::model::multi_polygon<polygon>;
  using mpoint = boost::geometry::model::multi_point<point>;
}

namespace spherical {
  using point = boost::geometry::model::d2::point_xy<double, boost::geometry::cs::spherical_equatorial<boost::geometry::degree> >;
  using linestring = boost::geometry::model::linestring<point>;
  using polygon = boost::geometry::model::polygon<point>;
  using multipolygon = boost::geometry::model::multi_polygon<polygon>;
  using mpoint = boost::geometry::model::multi_point<point>;
  const double earthRadiusMts = 6371000;
}

BOOST_AUTO_TEST_CASE(testUnionBug) {
  auto p1wkt = "POLYGON((-78.4072265625001 43.06652924482626,-78.4072265625 43.06740063068311,-78.4063141178686 43.06653210403569,-78.4072265625001 43.06652924482626))";
  auto p2wkt = "POLYGON((-78.55968743491499 43.06594969590624,-78.55036227331367 43.07380195109801,-78.53503704605811 43.08248347074284,-78.51769210872999 43.08880392487917,-78.49899564953199 43.09251971058174,-78.47966844278045 43.09348761253013,-78.46045580120891 43.09167037120638,-78.44209853911326 43.08713812460473,-78.42530412309867 43.08006566649393,-78.41071917768537 43.07072563376782,-78.40631359930124 43.06653210565861,-78.55968743491499 43.06594969590624))";

  {
    spherical::polygon p1,p2;
    boost::geometry::read_wkt(p1wkt, p1);
    boost::geometry::read_wkt(p2wkt, p2);

    spherical::multipolygon out;
    boost::geometry::union_(
          p1,
          p2,
          out);
    std::cout << "Spherical union: " << boost::geometry::wkt(out) << std::endl;
  }
  {
    cartesian::polygon p1,p2;
    boost::geometry::read_wkt(p1wkt, p1);
    boost::geometry::read_wkt(p2wkt, p2);

    cartesian::multipolygon out;
    boost::geometry::union_(
          p1,
          p2,
          out);
    std::cout << "Cartesian union: " << boost::geometry::wkt(out) << std::endl;
  }
}

when executed against Boost 1.66 the output is:

Spherical union: MULTIPOLYGON()
Cartesian union: MULTIPOLYGON(((-78.5597 43.0659,-78.5504 43.0738,-78.535 43.0825,-78.5177 43.0888,-78.499 43.0925,-78.4797 43.0935,-78.4605 43.0917,-78.4421 43.0871,-78.4253 43.0801,-78.4107 43.0707,-78.4063 43.0665,-78.5597 43.0659)))

Now I assume that union can only be empty if both constituting geometries are empty and under this assumption i take it this a bug. I will continue trying to diagnose and put more info, but as diagnostics are time consuming, I am submitting the issue as i see it now.

This the same general problem domain myself and @vosst have been working on (and earlier raised #470).

[remster]

I will just keep filling this comment up with my findings as i diagnose this.

I see that cartesian and spherical scenarios start to differ here where the cartesian finds no turns, whilst spherical does. Moreover turn_info_per_ring in the spherical case has, for both turns, has_blocked_turn set, which seems to influence ring selection here

[awulkiew]

Thanks for the info. I'll look into it next week.

[remster]

@awulkiew are there any news on this? I could have another go at diagnosing this, but the code is quite cryptic. Is it described somewhere? Does it implement a specific, well known algorithm? - so i can build a mental model of what it's trying to do.

[awulkiew]

@remster AFAIK the algorithm is Boost.Geometry-specific algorithm based on graph traversal. @barendgehrels was it taken from some paper or designed from scratch?

I can reproduce the bug and currently am working on the fix. I have a few possible fixes in mind (in spherical intersection strategy) but I want to debug more before I choose the best one because it's possible that the bug is located elsewhere (in overlay: handle_colocations, enrich_turns, traverse, etc).
Is this the only case where you get wrong result in spherical CS?
Is this something critical for you? How fast do you need the fix?

FYI, replacing the long segment in the second geometry with shorter ones seems to help, so e.g.:

boost::geometry::densify(p2, p3, 0.001); // for unit sphere
boost::geometry::union_(p1, p3, out); // instead of p2

[remster]

Thanks for the reply @awulkiew and for diving into this.

Is this the only case where you get wrong result in spherical CS?

Nope. Attaching more samples. I can't tell if they're down to the same/single defect or expose more. Note that we can only detect them in bulk by watching out for an empty union (which I assume - and given well formed inputs, is always an illegal output). We won't be able to do this with intersections.
out.txt

Is this something critical for you? How fast do you need the fix?

We're migrating our solution from gdal/geos to boost and we're doing this for two reasons:

1. header only vs a bulky library
2. one of our deliverables are mapbox tiles and we find it easier to generate them of off boost geometries - they're more transparent and flexible than gdal/geos. Doing so also looks about 2.2x quicker, which, to say the least, is important.

We've reached the functional parity and are discovering those bugs as we harden the solution. We obviously cannot complete the migration with these defects in - so we're having a stale branch that's getting gradually harder to keep in sync. If i were convinced that this is the only bug remaining, i could wait a month for a fix, but I am not convinced.

I would be happy to help diagnosing/fixing this - but I'd need a hand from you - to get me up to speed a bit.

[remster]

Found another possible problem.

namespace cartesian {
  using point = boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>;
  using linestring = boost::geometry::model::linestring<point>;
  using polygon = boost::geometry::model::polygon<point>;
  using multipolygon = boost::geometry::model::multi_polygon<polygon>;
  using mpoint = boost::geometry::model::multi_point<point>;
}

namespace spherical {
  using point = boost::geometry::model::d2::point_xy<double, boost::geometry::cs::spherical_equatorial<boost::geometry::degree> >;
  using linestring = boost::geometry::model::linestring<point>;
  using polygon = boost::geometry::model::polygon<point>;
  using multipolygon = boost::geometry::model::multi_polygon<polygon>;
  using mpoint = boost::geometry::model::multi_point<point>;
  const double earthRadiusMts = 6371000;
}

BOOST_AUTO_TEST_CASE(usaIntersectionTestCase) {
  auto p1wkt = "POLYGON((-180 -2.25,-180 80,4.5 80,4.5 -2.25,-180 -2.25))";
  auto p2wkt = "POLYGON((-125.0840938165823 48.18451089135078,-123.2484442951568 48.28402780487058,-123.3222396381058 49.00207166780857,-106.2294164816153 48.99936106028161,-94.95738881549684 49.37019439405368,-91.56749988130491 48.04377777876205,-88.36986182651864 48.30606296052222,-82.51861102256366 45.33863888434129,-82.67972207213995 41.67655556472526,-74.97255548445517 44.98347221572623,-71.50101930291031 45.01344579588383,-69.22446981284281 47.4598397141629,-67.79019834014535 47.06722478688557,-67.04184250949656 44.36199518537477,-70.23005912341451 43.2124078816444,-69.75267183371977 41.13341339615248,-73.71086681700677 40.38246941502378,-75.73925419455345 36.92336928786909,-75.32570262625269 35.10289000177789,-81.00038138644489 31.36986317745486,-80.21558018606714 24.87617100882401,-82.31642199133823 24.40399007010678,-81.37277484337314 25.35110770766313,-82.98991779576328 28.87424738341057,-86.63551567900944 30.19587268651689,-88.60142307258893 30.01998991835538,-89.40197910608261 28.705405956028,-94.4542152059084 29.30489831467678,-96.8117767445914 27.80056771982362,-97.42440989188675 25.84041996728249,-99.0845302237628 26.39797732094523,-101.4011198718317 29.77120000013631,-103.2862798604784 28.97691995351909,-106.4506798904261 31.76422996252971,-111.074824919583 31.33223902311842,-114.7199604399525 32.71865525370477,-117.4636899326359 32.58936795833381,-118.1537433047147 33.49758670948682,-118.5200689113933 32.61376470726387,-119.9333407371885 33.34100477748434,-118.6596572250774 33.83057945961541,-120.7240783642898 33.9956452157932,-124.7155563780111 40.37995067155807,-125.0840938165823 48.18451089135078))";

  {
    spherical::polygon p1,p2;
    boost::geometry::read_wkt(p1wkt, p1);
    boost::geometry::read_wkt(p2wkt, p2);

    spherical::multipolygon out;
    boost::geometry::intersection(p1,p2,out);
    std::cout << "Spherical intersection: " << boost::geometry::wkt(out) << std::endl;
  }
  {
    cartesian::polygon p1,p2;
    boost::geometry::read_wkt(p1wkt, p1);
    boost::geometry::read_wkt(p2wkt, p2);

    cartesian::multipolygon out;
    boost::geometry::intersection(p1,p2,out);
    std::cout << "Cartesian intersection: " << boost::geometry::wkt(out) << std::endl;
  }
}

prints

Spherical intersection: MULTIPOLYGON()
Cartesian intersection: MULTIPOLYGON(((-125.084 48.1845,-123.248 48.284,-123.322 49.0021,-106.229 48.9994,-94.9574 49.3702,-91.5675 48.0438,-88.3699 48.3061,-82.5186 45.3386,-82.6797 41.6766,-74.9726 44.9835,-71.501 45.0134,-69.2245 47.4598,-67.7902 47.0672,-67.0418 44.362,-70.2301 43.2124,-69.7527 41.1334,-73.7109 40.3825,-75.7393 36.9234,-75.3257 35.1029,-81.0004 31.3699,-80.2156 24.8762,-82.3164 24.404,-81.3728 25.3511,-82.9899 28.8742,-86.6355 30.1959,-88.6014 30.02,-89.402 28.7054,-94.4542 29.3049,-96.8118 27.8006,-97.4244 25.8404,-99.0845 26.398,-101.401 29.7712,-103.286 28.9769,-106.451 31.7642,-111.075 31.3322,-114.72 32.7187,-117.464 32.5894,-118.154 33.4976,-118.52 32.6138,-119.933 33.341,-118.66 33.8306,-120.724 33.9956,-124.716 40.38,-125.084 48.1845)))

It's this case:

[awulkiew]

Thanks. I'll look into this as well. It looks like another issue. The result should be the US.

But that's interesting. Just to be sure we look at this problem the same way. In spherical this is not how the tile's Polygon (p1) looks like. Edges of Polygons are defined by great circles in Boost.Geometry so horizontal edges of this tile are not parallels. So e.g. the intersection point of two segments will be different in spherical and cartesian because the line interpolation is different. Moreover because the edges of Polygons are defined by great circles it is not possible to represent the intersection of this tile's shape (something which in Boost.Geometry could be considered as Box with 2 edges being parallels) and a Polygon. And before you try to represent the tile as Box in the code. This will probably compile but currently a Box is internally treated as a Polygon in all coordinate systems which is wrong in non-cartesian ones.

densify() could be used on this tile Polygon to see what it looks like on a sphere.

[remster]

Sure.. the vertical edges of this tile are of course non-parallel in geometric sense, but they run along their respective meridians so i can't imagine the picture changing post-densify and when projected on a web mercator.
densify appears on trunk and I am working of off 1.66 so can't reach it.

I assume we agree that the intersection should not be empty in the spherical coordinate system.

[awulkiew]

I'm talking about horizontal edges. When represented in Polygon they are not going along parallels so lines of the same latitude (https://en.wikipedia.org/wiki/Parallel_(latitude)). Instead they are going along great circles (https://en.wikipedia.org/wiki/Great_circle). One exception is equator.

Vertical edges are going along meridians which are equivalent to great circles.

Actually I'm taking back my previous statement "The result should be the US.". The Polygon is: POLYGON((-180 -2.25,-180 80,4.5 80,4.5 -2.25,-180 -2.25)). Consider one of the horizontal edges: -180 80,4.5 80. This edge will be defined by great circle crossing both of these points. Now you have 2 possibilities going west or east. The shorter path between points is choosen by Boost.Geometry. So this segment goes from longitude -180 to longitude 4.5 on the east hemisphere (wrt. prime meridian). So from east Asia to the west, through the whole Asia to west Europe. So this Polygon does not cover the US.

And taking into account all of the above this Polygon is invalid because it's CCW, not CW (it has to be reversed). Btw, Boost.Geometry considers longitude 180 as valid, not -180. -180 will work but all operations will generate 180 so if you don't want to have them mixed use 180.

This is how this Polygon looks like before (red) and after (green) densification in a cartesian-like spherical coordinate space:

[remster]

Very interesting - i guess i assumed the winding order will be a hint stronger than edge lengths.
So all i have to do (if i understand correctly) is to densify the big tile poly to have its edges sufficiently shorter (certainly less than the half distance of the great circle). I think i gonna have to start working from trunk, as, correct me if i am wrong - densify is not released.

Let's put that last issue aside for now - i'll test at work tomorrow. I still hope you can help us with the empty unions.

BTW, does a CCW polygon become invalid (it has no defined meaning) or is it interpreted as a hole in the sphere (i.e.: it occupies the entirety of the coordinate system minus the CCW ring)? I certainly saw evidence to suggest this but it could be a side-effect of the poly's invalidness-ness.

You said denisfy goes around the original problem i reported, I will scratch my head on that, but i guess it's unrelated as these geometries are small.

[awulkiew]

So all i have to do (if i understand correctly) is to densify the big tile poly to have its edges sufficiently shorter

Note that if you want to represent this tile then you don't have to use densify(). You can simply generate more points on the same parallel (the same latitude) on the horizontal edges. The more points you create the more accurate the tile is represented. However it also means that more points will be in the output of set operation for this Polygon representing the tile and some other Polygon intersecting the edge of the tile. In this case one Point is needed on each horizontal edge at minimum:

I think i gonna have to start working from trunk, as, correct me if i am wrong - densify is not released.

Yes, it is a new feature in 1.67.

I still hope you can help us with the empty unions.

Yes, I'm working on it.

does a CCW polygon become invalid (it has no defined meaning) or is it interpreted as a hole in the sphere (i.e.: it occupies the entirety of the coordinate system minus the CCW ring)?

It depends on algorithm. Some algorithms works with such Polygons but in general they are considered as invalid at the moment. As you mentioned they are used to represent holes. It's because the same mechanisms are used for cartesian and non-cartesian Rings when some algorithm has to calculate the order. If both orders were correct then it would be impossible to calculate one of the two possible orders while analysing only the data, so Points in a Ring. This is the case in set operations for example.

You said denisfy goes around the original problem i reported, I will scratch my head on that, but i guess it's unrelated as these geometries are small.

The original problem is caused by numerical issues inside the algorithm. At some point the algorithm has to check at which place on a segment the intersection point lies. If the lengths of segments are different then calculations are performed for different numbers.