diff options
| -rw-r--r-- | lib/great_circle_arc.py | 23 | ||||
| -rw-r--r-- | lib/polygon.py | 26 |
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) |