Geometry Boost List: Re: [geometry] Spherical area

Subject: Re: [geometry] Spherical area
From: Barend Gehrels (barend_at_[hidden])
Date: 2016-08-10 05:35:17

Hi Vissarion,

Thanks and sorry for the long delay.

Op 7-7-2016 om 20:20 schreef Vissarion Fysikopoulos:
> Hi all,
>
> I am working on the computation of the spherical area of a polygon. I 
> would like to discuss two topics.
>
>
> A. Currently Boost.Geometry has a strategy called area_huiller (there 
> is a typo, the right name is Huilier). Basically, this strategy adds 
> the signed areas of the triangles formed by the North pole and a 
> segment of the polygon.

It is OK to correct the name.

>
> Apart from the typo there are two more serious problems about this 
> strategy.
>
> 1) Accuracy problems because of badly conditioned triangles (cf. 
> https://svn.boost.org/trac/boost/ticket/11930)
>
> 2) Incorrect results in corner cases: polygons that contain a pole, 
> polygons that has a pole as a vertex (cf. 
> https://svn.boost.org/trac/boost/ticket/11931)
>
> Problem 1) is fixed by using another strategy (let me call it 
> area_trapezoidal) that adds the areas of the "spherical trapezoids" 
> that formed by each segment of the polygon tow meridian segments 
> starting from segment's endpoints and end at the equator, and the 
> corresponding equator segment. See the last formula in 
> https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess
>
> Problem 2) can be solved independently in both strategies, by keeping 
> track whether the polygon encircles the South pole and then subtract a 
> quantity---4*\pi or 2*\pi depending on the strategy---from the result.
>
> For the case of polygons that have South pole as a vertex (sub-case of 
> Problem 2), area_huiller strategy should do more work (either switch 
> the computation using the North pole or count the angle of the polygon 
> on South pole and subtract the corresponding quantity).
>
> The first question here is: does it worth the effort to implement 
> fixes for area_huiller since the trapezoidal one seems to be supreme? 
> Do you have any test case where you need a formula such as 
> area_huiller, i.e. it assumes that we decompose the polygon to 
> triangles? My opinion is the area_huiller is redundant and we can 
> remove it.

I'm OK with replacing it with area_trapezoidal too, as long as it is as 
fast. If it is considerably slower, I prefer to keep both to give 
alternatives for fast computation (but with unsolved pole behaviour) and 
for correct computation for all cases, including geographic.

>
> Another advantage of area_trapezoidal is that it can be used as a 
> basis for geographic area, which will be my next step.
>
>
> B. There is a natural question on the semantics of the spherical 
> polygon area. Given a polygon on a sphere is not clear if we refer to 
> the polygon with the smallest enclosing area or it's complement.

It is indeed ambiguous. See also below.

> Also, I guess that we assume that the segments of a polygon are 
> geodesics to avoid cases of polygons with segments that wrap around 
> the globe. For example POLYGON((1 1, -1 1, -1 -1, 1 -1, 1 1)) can be 
> assumed to wrap around the sphere or not. With the assumption of 
> geodesic segments this is uniquely defined.

I am OK with that.

>
> This issue is connected to whether we allow a polygon to have negative 
> area. A polygon admits two orientations. The one corresponds to 
> positive area. The other either corresponds to the same value of area 
> negated or to the area of the complement of the polygon. I would go 
> for the first definition since it seems that will be more compatible 
> to cartesian case.

It is also related to interior rings. So a polygon is (depending on the 
template) clockwise oriented and then has a positive area. If it is 
inverted, its area is negative but for a polygon that should not exist. 
For a ring though, all interior rings should have a negative area.
So also on the sphere this should have the same definition. But there 
indeed there is another interpretation possible, as you describe. We 
could allow that, because on cartesian a negative polygon "should not 
exist" as I mentioned, but on the sphere it could exist and be the 
complement, still having a positive area but spanning more than a 
hemisphere.

>
> From your experience with Boost.Geometry what do you think is the best 
> semantics?

Regards, Barend