From b6749b81e09abad7944a4af21417e23cd8ad3460 Mon Sep 17 00:00:00 2001 From: mdroe Date: Mon, 4 Jun 2012 18:58:11 +0000 Subject: 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 --- lib/great_circle_arc.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) (limited to 'lib/great_circle_arc.py') 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 -- cgit