spacetelescope / spherical_geometry

A Python package for handling spherical polygons that represent arbitrary regions of the sky
http://spherical-geometry.readthedocs.io/
62 stars 31 forks source link

Assertion error when computing intersection of nearly overlapping polygons #168

Open mcara opened 5 years ago

mcara commented 5 years ago

When computing intersections of almost overlapping polygons. Example:

>>> import numpy as np
>>> from spherical_geometry.polygon import SphericalPolygon
>>> p1 = SphericalPolygon(
...     np.array([[-0.1094946215827374, -0.8592766830993238, -0.499654390280199 ],
...               [-0.1089683641318892, -0.8595220381654031, -0.4993473355555343],
...               [-0.108610535224965 , -0.8593183788298407, -0.4997756250993051],
...               [-0.1091500557209236, -0.8590667764452905, -0.5000905307482003],
...               [-0.1094946215827374, -0.8592766830993238, -0.499654390280199 ]]),
...     np.array([-0.1090595793730483, -0.8592979843505629, -0.4997128998115153])
... )
>>> p2 = SphericalPolygon(
...     np.array([[-0.1094946213367254, -0.8592766831114167, -0.4996543903133135],
...               [-0.1089683641834766, -0.859522038038747 , -0.4993473357622887],
...               [-0.1086105354789061, -0.8593183788183577, -0.4997756250638628],
...               [-0.109150055669766 , -0.8590667765760884, -0.5000905305346783],
...               [-0.1094946213367254, -0.8592766831114167, -0.4996543903133135]]),
...     np.array([-0.1090595793730483, -0.8592979843505629, -0.4997128998115153])
... )
>>> p1.intersection(p2)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/.../lib/python3.7/site-packages/spherical_geometry/polygon.py", line 1198, in intersection
    subpolygons = polya.intersection(polyb)
  File "/.../lib/python3.7/site-packages/spherical_geometry/polygon.py", line 657, in intersection
    return g.intersection()
  File "/.../lib/python3.7/site-packages/spherical_geometry/graph.py", line 437, in intersection
    self._sanity_check("intersection - find all intersections")
  File "/.../lib/python3.7/site-packages/spherical_geometry/graph.py", line 390, in _sanity_check
    assert len(node._edges) >= 2
AssertionError
nihargupte-ph commented 4 years ago

I am experiencing the same error, a potential temporary workaround, if you are only concerned with the area, for example, is projecting the spherical polygon to a shapely object and find the intersection that way.

``def spherical_poly_to_poly(poly):
    """ Converts a spherical polygon to a lon lat polygon """
    radec = list(poly.to_radec())
    lons, lats = radec[0][0], radec[0][1]
    poly = geometry.Polygon(list(zip(lons, lats)))
    return poly

I think part of the problem is some sort of self intersection in the polygons? When I convert to shapely objects I get errors even through shapely intersections but only for the ones that have assertion errors SphericalPolyon. A fix for this is shapely is to simply buffer the polygon --> p = p.buffer(0). I'll let you know if I find a fix