Use an entirely different method for area computation. It first transforms to a 2D plane by interpolating the great circle arcs through a Lambert azimuthal equal area projection. Then a standard 2D polygon area method is used.

Make sure the polygons always go clockwise to aid in area computation.

Fix union calculation -- it should be removing interior edges, not interior nodes.

Fix some numerical out-of-range problems in great_circle_arc.py

Remove the serial method for multi_union -- it no longer generates polygons that are compatible with calculating the area.

Make debugging images more explanatory by putting a meaningful title at the top.


git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci_python/branches/sphere@17200 fe389314-cf27-0410-b35b-8c050e845b92

Former-commit-id: fb841a481025150f631be6d64cf995dda592ecbc

1 files changed, 6 insertions, 2 deletions

diff --git a/lib/great_circle_arc.py b/lib/great_circle_arc.py
index 682f9d9..f70d16e 100644
--- a/lib/great_circle_arc.py
+++ b/lib/great_circle_arc.py
@@ -218,6 +218,7 @@ def length(A, B, degrees=True):
     B = np.asanyarray(B)
 
     dot = inner1d(A, B)
+    dot = np.clip(dot, -1.0, 1.0)
     with np.errstate(invalid='ignore'):
         result = np.arccos(dot)
 
@@ -350,7 +351,10 @@ def interpolate(A, B, steps=50):
     t = np.linspace(0.0, 1.0, steps, endpoint=True).reshape((steps, 1))
 
     omega = length(A, B, degrees=False)
-    sin_omega = np.sin(omega)
-    offsets = np.sin(t * omega) / sin_omega
+    if omega == 0.0:
+        offsets = t
+    else:
+        sin_omega = np.sin(omega)
+        offsets = np.sin(t * omega) / sin_omega
 
     return offsets[::-1] * A + offsets * B