polygon.SphericalPolygon.convex_hull() sometimes fails with a "ValueError: Null vector"

[pforshay]

While generating spatial footprints involving multiple polygons, we've run into sporadic cases where trying to generate the convex hull results in a "ValueError: Null vector." Our processing is set up to close all generated polygons by repeating the first coordinate pair again at the end. While investigating one particular case where the closed polygons resulted in this ValueError, we found that removing the repeated coordinates resulted in a successful convex hull calculation. However, we do not routinely run into these errors with other closed multi-polygon shapes, so it seems to be something about these specific cases that causes trouble for the convex hull.

Below is a bit of code that should replicate the error we are seeing, including the steps we take to get from a set of ra/dec coordinates to the convex hull calculation, along with the full output from the included fail case:

import numpy as np
from spherical_geometry import polygon, vector

def get_convex_hull(vx, vy):
    x, y, z = vector.lonlat_to_vector(vx, vy)
    points = np.dstack((x, y, z))[0]
    return polygon.SphericalPolygon.convex_hull(points)

def fail_case():
    # Two polygons with repeated closing coords, which fails in convex_hull
    vx = [175.580350, 175.616304, 175.636560, 175.600604, 175.580350, 175.600881, 175.636915, 175.656911, 175.620856, 175.600881]
    vy = [30.338842, 30.372084, 30.357442, 30.324553, 30.338842, 30.324353, 30.357186, 30.342694, 30.310229, 30.324353]
    return vx, vy

def pass_case_open():
    # Same two polygons as fail_case() without the closing coords
    vx = [175.580350, 175.616304, 175.636560, 175.600604, 175.600881, 175.636915, 175.656911, 175.620856]
    vy = [30.338842, 30.372084, 30.357442, 30.324553, 30.324353, 30.357186, 30.342694, 30.310229]
    return vx, vy

def pass_case_closed():
    # Two polygons with repeated closing coords, which work with convex_hull
    vx = [219.475836, 219.43058, 219.431779, 219.476735, 219.475836, 219.476753, 219.431842, 219.433006, 219.477659, 219.476753]
    vy = [-1.725932, -1.721104, -1.698335, -1.70337, -1.725932, -1.703015, -1.69796, -1.675434, -1.680744, -1.703015]
    return vx, vy

if __name__ == '__main__':
    cases = [pass_case_open, pass_case_closed, fail_case]
    hull = None
    for c in cases:
        try:
            hull = get_convex_hull(*c())
            print(f"--- {c.__name__} succeeded:")
        except ValueError as err:
            hull = err
            print(f"--- {c.__name__} failed:")
        print(hull)

Error output:

File /site-packages/spherical_geometry/polygon.py:984, in SphericalPolygon.convex_hull(cls, points)
    970 @classmethod
    971 def convex_hull(cls, points):
    972     r"""
    973     Create a new SphericalPolygon from the convex hull of a
    974     list of points.
   (...)
    982     polygon : SphericalPolygon object
    983     """
--> 984     polygon = SingleSphericalPolygon.convex_hull(points)
    985     return cls((polygon,))

File /site-packages/spherical_geometry/polygon.py:350, in SingleSphericalPolygon.convex_hull(cls, points)
    347 # Sort points around extreme point by angle from true north
    349 north = [0., 0., 1.]
--> 350 ang = great_circle_arc.angle(north, extreme, points)
    351 pt = [points[i,:] for i in range(points.shape[0])]
    353 duo = list(zip(pt, ang))

File /site-packages/spherical_geometry/great_circle_arc.py:322, in angle(A, B, C, degrees)
    298 """
    299 Returns the angle at *B* between *AB* and *BC*.
    300 
   (...)
    319    polygon.  Graphics Gems IV.  1994.  Academic Press.
    320 """
    321 if HAS_C_UFUNCS:
--> 322     angle = math_util.angle(A, B, C)
    323 else:
    324     A = np.asanyarray(A)

ValueError: Null vector.

[pforshay]

Adding package versions:

numpy                     1.24.3                   pypi_0    pypi
python                    3.11.0          he550d4f_1_cpython    conda-forge
spherical-geometry        1.2.23                   pypi_0    pypi

[mcara]

@pforshay Thanks for reporting this.