Negative Polygon Area

[sfinkens]

spherical_geometry.get_polygon_area() computes negative values near the poles. The result also depends on the order of the given corners:

from pyresample.spherical_geometry import Coordinate, get_polygon_area

ll = Coordinate(0, 89.925)
ul = Coordinate(0, 89.975)
ur = Coordinate(0.05, 89.975)
lr = Coordinate(0.05, 89.925)


print get_polygon_area([ll, ul, ur, lr])  # -0.00043633152367
print get_polygon_area([ll, ul, lr, ur])  # 9.96988713808e-10

Is there an expected order of the corners?

[mraspaud]

Yes, the order of the vertices is important. A polygon is defined with vertices in a clockwise order around it IIRC.

But the area should never be negative, thanks for reporting this.

[sfinkens]

Thanks for the clarification! According to this source, the area of a spherical triangle with angles A,B,C is

R^2*((A+B+C) - PI)

However, the current implementation in get_polygon_area() is

area += R**2 * e__ - math.pi

which is equivalent to

R^2*(A+B+C) - PI

Maybe that is a problem?

[sfinkens]

Ah no, as long as R=1 there is no difference. But for the future it might be good to fix that:

def get_polygon_area(corners):
    """Get the area of the convex area defined by *corners*.
    """
    # We assume the earth is spherical !!!
    # Should be the radius of the earth at the observed position
    R = 1

    c1_ = corners[0]
    area = 0

    for idx in range(1, len(corners) - 1):
        b1_ = Arc(c1_, corners[idx])
        b2_ = Arc(c1_, corners[idx + 1])
        b3_ = Arc(corners[idx], corners[idx + 1])
        e__ = (abs(b1_.angle(b2_)) +
               abs(b2_.angle(b3_)) +
               abs(b3_.angle(b1_)))
        area += e__ - math.pi
    return R**2 * area

Shall I send a pull request?

[mraspaud]

A PR would be great.

Thanks !

[sfinkens]

I think I also found the reason for negative polygon areas: In spherical_geometry.Arc.angle() the angle between two arcs is computed using the dot product of normal vectors ua_ and ub_:

val = ua_.dot(ub_) / (ua_.norm() * ub_.norm())
if abs(val - 1) < EPSILON:
    angle = 0
elif abs(val + 1) < EPSILON:
    angle = math.pi
else:
    angle = math.acos(val)

If the angle is very small, i.e. the dot product is almost 1, the angle is reset to zero in the subsequent if block. That causes the sum e__ of angles in the spherical triangle to be slightly less than pi which finally leads to a negative polygon area.

I assume the EPSILON criteria were implemented to prevent dot products like 1.00...1 from triggering math domain errors in math.acos()? In that case I would propose the following solution which only takes dot products > 1 and < -1 into account:

if 0 <= val - 1 < EPSILON:
    angle = 0
elif -EPSILON < val + 1 <= 0:
    angle = math.pi
else:
    angle = math.acos(val)

[mraspaud]

Well actually everything outside of [-1, 1] would trigger an error, so maybe we could simply write

if val > 1: angle = 0
etc... ?

[sfinkens]

That's true. But it would mask a possible bug in the dot product computation: Even if the dot product yielded a completely unrealistic value, the angles would still be valid.

[mraspaud]

Good point.

Also, I remember now the use of EPSILON: it's also to be able to detect colinearity, so smaller angles are snapped to 0...

[sfinkens]

I see. What about adding a snap argument to the angle() method,

if snap:
    if abs(val - 1) < EPSILON:
        angle = 0
    elif abs(val + 1) < EPSILON:
        angle = math.pi
    else:
        angle = math.acos(val)
else:
    if 0 <= val - 1 < EPSILON:
        angle = 0
    elif -EPSILON < val + 1 <= 0:
        angle = math.pi
    else:
        angle = math.acos(val)

and calling it with snap=False in get_polygon_area()?

[mraspaud]

Sounds good 👍

[sfinkens]

The last two tests on overlap rate would then fail:

======================================================================
FAIL: Test how much two areas overlap.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/sfinkens/code/pyresample/pyresample/test/test_spherical_geometry.py", line 154, in test_overlap_rate
    self.assertAlmostEqual(area1.overlap_rate(area2), 0.5, 2)
AssertionError: 0.5086952608806948 != 0.5 within 2 places

======================================================================
FAIL: Test how much two areas overlap.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/sfinkens/code/pyresample/pyresample/test/test_spherical_geometry.py", line 155, in test_overlap_rate
    self.assertAlmostEqual(area2.overlap_rate(area1), 0.068, 3)
AssertionError: 0.06850857188651323 != 0.068 within 3 places

Although the computed values are very close. Original values are 0.499837070923 and 0.0679026012212. Rounding to (1,2) instead of (2,3) digits would resolve the issue, but I don't know whether this is correct.

[mraspaud]

Looking at the test code, I don't think there is a reason 0.499837070923 should be more correct than 0.5086952608806948. And since you removed the snapping, I think the latter should be preferred.
What about changing the test values to 0.509 and 0.0685 ?

[sfinkens]

Ok!

[mraspaud]

Merged, thanks for bringing this up !

[sfinkens]

Happy to help. Thanks for the great work you have done already!