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, 31 insertions, 3 deletions

diff --git a/lib/vector.py b/lib/vector.py
index 7323c92..5d26161 100644
--- a/lib/vector.py
+++ b/lib/vector.py
@@ -198,9 +198,6 @@ def rotate_around(x, y, z, u, v, w, theta, degrees=True):
     if degrees:
         theta = np.deg2rad(theta)
 
-    u2 = u**2
-    v2 = v**2
-    w2 = w**2
     costheta = np.cos(theta)
     sintheta = np.sin(theta)
     icostheta = 1.0 - costheta
@@ -211,3 +208,34 @@ def rotate_around(x, y, z, u, v, w, theta, degrees=True):
     Z = (-w*det)*icostheta + z*costheta + (-v*x + u*y)*sintheta
 
     return X, Y, Z
+
+
+def equal_area_proj(x, y, z):
+    """
+    Transform the coordinates to a 2-dimensional plane using the
+    Lambert azimuthal equal-area projection.
+
+    Parameters
+    ----------
+    x, y, z : scalars or 1-D arrays
+        The input vectors
+
+    Returns
+    -------
+    X, Y : tuple of scalars or arrays of the same length
+
+    Notes
+    -----
+
+    .. math::
+
+        X = \sqrt{\frac{2}{1-z}}x
+
+    .. math::
+
+        Y = \sqrt{\frac{2}{1-z}}y
+    """
+    scale = np.sqrt(2.0 / (1.0 - z))
+    X = scale * x
+    Y = scale * y
+    return X, Y