SphericalGeometry overlap fails

[mcdigman]

SphericalGeometry overlap, intersection, and union all fail with an AssertionError when given nested geometries with 1 identical edge. The following code demonstrates the issue:

from spherical_geometry.polygon import SphericalPolygon
import spherical_geometry.vector as sgv
import numpy as np

#set up nested polygons with an overlapping edge
bounding_theta1 = np.array([-np.pi/4.,0.,0.0,-np.pi/4.,-np.pi/4.])
bounding_theta2 = np.array([-np.pi/4.,0.,0.0,-np.pi/4.,-np.pi/4.])
bounding_phi1 = np.array([0.,0.,np.pi/4.,np.pi/4.,0.])
bounding_phi2 = np.array([0.,0.,np.pi/2.,np.pi/2.,0.])
theta_in = -np.pi/8.
phi_in = np.pi/8.

bounding_xyz1 = np.array(sgv.radec_to_vector(bounding_phi1,bounding_theta1,degrees=False)).T
bounding_xyz2 = np.array(sgv.radec_to_vector(bounding_phi2,bounding_theta2,degrees=False)).T
inside_xyz = np.array(sgv.radec_to_vector(phi_in,theta_in,degrees=False))

sp_poly1 = SphericalPolygon(bounding_xyz1,inside=inside_xyz)
sp_poly2 = SphericalPolygon(bounding_xyz2,inside=inside_xyz)

#should give 1., gives error
print sp_poly1.overlap(sp_poly2)

[bernie-simon]

I'm not getting the error, instead I get the answer 0.0. But there are two problems here. First, since the one polygon is entirely within the second the answer should not be zero. Looking at the source code, SphericalPolygon seems to handle the case where one polygon intersects another, but not the case where one polygon is entirely contained in another. I think there needs to be additional ocde to check for the contained case and handle it properly. Second, although I did not get the error you report, I strongly believe your problem is caused by rounding error. The code that computes the intersection between two lines is most inaccurate when the two lines are coincident or nearly coincident, because the cross product is least accurate in this case. I've thought about modifying the intersection code to test for the case of coincident lines and use different code for that case.

It would help if you sent me the full printout from the assertion error that you get, as I cannot duplicate the problem myself.

[mcdigman]

It produces this error and the attached plot:

Traceback (most recent call last): File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/spherical_geometry/graph.py", line 392, in _sanity_check assert len(node._edges) >= 2 AssertionError

I can also replicate the behavior of failing and return 0.0 without error, for example by modifying this line in the original example:
bounding_theta1 = np.array([-np.pi/4.,0.,0.0+0.001,-np.pi/4.,-np.pi/4.])

[mcara]

With the latest release I am getting the expected value 1. I think this is no longer an issue.