Fix drawing -- basemap can not handle drawing great circle arcs that go around the edges of the domain. The solution is to do the interpolation in the quaternion space and then pass little snippets of that to basemap.

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

Former-commit-id: dc807419ae1618c42af6434f6e63c7f0cb3bf36e

2 files changed, 38 insertions, 11 deletions

diff --git a/lib/great_circle_arc.py b/lib/great_circle_arc.py
index 056857d..93d9ad2 100644
--- a/lib/great_circle_arc.py
+++ b/lib/great_circle_arc.py
@@ -315,3 +315,26 @@ def midpoint(A, B):
     l = np.sqrt(np.sum(P * P, axis=-1))
     l = np.expand_dims(l, 2)
     return P / l
+
+
+def interpolate(A, B, steps=50):
+    """
+    Interpolate along the great circle arc.
+
+    Parameters
+    ----------
+    A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
+        The endpoints of the great circle arc.  It is assumed that
+        these points are already normalized.
+
+    steps : int
+        The number of interpolation steps
+
+    Returns
+    -------
+    array : (*x*, *y*, *z*) triples
+        The points interpolated along the great circle arc
+    """
+    t = np.linspace(0, 1.0, steps, endpoint=True).reshape((steps, 1))
+
+    return (t * A) + ((1.0 - t) * B)
diff --git a/lib/polygon.py b/lib/polygon.py
index 79f1d64..6241416 100644
--- a/lib/polygon.py
+++ b/lib/polygon.py
@@ -53,12 +53,11 @@ __all__ = ['SphericalPolygon']
 
 class SphericalPolygon(object):
     ur"""
-    Polygons are represented by both a set of points (in
-    Cartesian (*x*, *y*, *z*) normalized on the unit sphere),
-    and an inside point.  The inside point is necessary, because
-    both the inside and outside of the polygon are finite areas
-    on the great sphere, and therefore we need a way of
-    specifying which is which.
+    Polygons are represented by both a set of points (in Cartesian
+    (*x*, *y*, *z*) normalized on the unit sphere), and an inside
+    point.  The inside point is necessary, because both the inside and
+    outside of the polygon are finite areas on the great sphere, and
+    therefore we need a way of specifying which is which.
     """
 
     def __init__(self, points, inside):
@@ -668,11 +667,16 @@ class SphericalPolygon(object):
         if not len(plot_args):
             plot_args = {'color': 'blue'}
         points = self._points
-        ra, dec = vector.vector_to_radec(
-            points[:, 0], points[:, 1], points[:, 2],
-            degrees=True)
-        for r0, d0, r1, d1 in zip(ra[0:-1], dec[0:-1], ra[1:], dec[1:]):
-            m.drawgreatcircle(r0, d0, r1, d1, **plot_args)
+
+        for A, B in zip(points[0:-1], points[1:]):
+            length = great_circle_arc.length(A, B, degrees=True)
+            interpolated = great_circle_arc.interpolate(A, B, length * 4)
+            ra, dec = vector.vector_to_radec(
+                interpolated[:, 0], interpolated[:, 1], interpolated[:, 2],
+                degrees=True)
+            for r0, d0, r1, d1 in zip(ra[0:-1], dec[0:-1], ra[1:], dec[1:]):
+                m.drawgreatcircle(r0, d0, r1, d1, **plot_args)
+
         ra, dec = vector.vector_to_radec(
             *self._inside, degrees=True)
         x, y = m(ra, dec)