improved robustness of point in polygon-on-sphere tests

doc/source/library/reference/orthodrome.rst

@@ -3,7 +3,6 @@
 
 .. automodule:: pyrocko.orthodrome
     :members:
-    :undoc-members:
 
 
 

src/orthodrome.py

@@ -1200,50 +1200,107 @@ def stereographic_poly(points):
     return points_out
 
 
-def contains_point(polygon, point):
-    rot = rot_to_00(point[0], point[1])
-    points_xyz = latlon_to_xyz(polygon)
-    points_rot = num.dot(rot, points_xyz.T).T
-    groups = spoly_cut([points_rot], axis=0)
-    for group in groups:
-        for points_g in group:
-            try:
-                points2 = stereographic_poly(points_g)
-                p = Path(points2, closed=True)
-                if p.contains_point((0., 0.)):
-                    return True
+def gnomonic_x(points, cutoff=0.01):
+    points_out = points[:, 1:].copy()
+    if num.any(points[:, 0] < cutoff):
+        raise Farside()
+
+    factor = 1.0 / points[:, 0]
+    points_out *= factor[:, num.newaxis]
+    return points_out
 
-            except Farside:
-                pass
 
-    return False
+def cneg(i, x):
+    if i == 1:
+        return x
+    else:
+        return num.logical_not(x)
 
 
 def contains_points(polygon, points):
+    '''
+    Test which points are inside polygon on a sphere.
+
+    :param polygon: Point coordinates defining the polygon [deg].
+    :type polygon: :py:class:`numpy.ndarray` of shape (N, 2), second index
+        0=lat, 1=lon
+    :param points: Coordinates of points to test [deg].
+    :type points: :py:class:`numpy.ndarray` of shape (N, 2), second index
+        0=lat, 1=lon
+
+    :returns: Boolean mask array.
+    :rtype: :py:class:`numpy.ndarray` of shape (N,)
+
+    The inside of the polygon is defined as the area which is to the left hand
+    side of an observer walking the polygon line, points in order, on the
+    sphere. Lines between the polygon points are treated as great circle paths.
+    The polygon may be arbitrarily complex, as long as it does not have any
+    crossings or thin parts with zero width. The polygon may contain the poles
+    and is allowed to wrap around the sphere multiple times.
+
+    The algorithm works by consecutive cutting of the polygon into (almost)
+    hemispheres and subsequent Gnomonic projections to perform the
+    point-in-polygon tests on a 2D plane.
+    '''
+
+    and_ = num.logical_and
+
     points_xyz = latlon_to_xyz(points)
-    center_xyz = num.mean(points_xyz, axis=0)
+    mask_x = 0. <= points_xyz[:, 0]
+    mask_y = 0. <= points_xyz[:, 1]
+    mask_z = 0. <= points_xyz[:, 2]
 
-    assert num.all(
-        distances3d(points_xyz, center_xyz[num.newaxis, :]) < 1.0)
+    result = num.zeros(points.shape[0], dtype=num.int)
 
-    lat, lon = xyz_to_latlon(center_xyz)
-    rot = rot_to_00(lat, lon)
+    for ix in [-1, 1]:
+        for iy in [-1, 1]:
+            for iz in [-1, 1]:
+                mask = and_(
+                    and_(cneg(ix, mask_x), cneg(iy, mask_y)),
+                    cneg(iz, mask_z))
 
-    points_rot_xyz = num.dot(rot, points_xyz.T).T
-    points_rot_pro = stereographic(points_rot_xyz)
+                center_xyz = num.array([ix, iy, iz], dtype=num.float)
 
-    poly_xyz = latlon_to_xyz(polygon)
-    poly_rot_xyz = num.dot(rot, poly_xyz.T).T
-    groups = spoly_cut([poly_rot_xyz], axis=0)
-    result = num.zeros(points.shape[0], dtype=num.int)
-    for group in groups:
-        for poly_rot_group_xyz in group:
-            try:
-                poly_rot_group_pro = stereographic_poly(poly_rot_group_xyz)
-                result += path_contains_points(
-                    poly_rot_group_pro, points_rot_pro)
+                lat, lon = xyz_to_latlon(center_xyz)
+                rot = rot_to_00(lat, lon)
 
-            except Farside:
-                pass
+                points_rot_xyz = num.dot(rot, points_xyz[mask, :].T).T
+                points_rot_pro = gnomonic_x(points_rot_xyz)
+
+                offset = 0.01
+
+                poly_xyz = latlon_to_xyz(polygon)
+                poly_rot_xyz = num.dot(rot, poly_xyz.T).T
+                poly_rot_xyz[:, 0] -= offset
+                groups = spoly_cut([poly_rot_xyz], axis=0)
+
+                for poly_rot_group_xyz in groups[1]:
+                    poly_rot_group_xyz[:, 0] += offset
+
+                    poly_rot_group_pro = gnomonic_x(
+                        poly_rot_group_xyz)
+
+                    if circulation(poly_rot_group_pro) > 0:
+                        result[mask] += path_contains_points(
+                            poly_rot_group_pro, points_rot_pro)
+                    else:
+                        result[mask] -= path_contains_points(
+                            poly_rot_group_pro, points_rot_pro)
 
     return result.astype(num.bool)
+
+
+def contains_point(polygon, point):
+    '''
+    Test if point is inside polygon on a sphere.
+
+    :param polygon: Point coordinates defining the polygon [deg].
+    :type polygon: :py:class:`numpy.ndarray` of shape (N, 2), second index
+        0=lat, 1=lon
+    :param point: Coordinates ``(lat, lon)`` of point to test [deg].
+
+    Convenience wrapper to :py:func:`contains_points` to test a single point.
+    '''
+
+    return bool(
+        contains_points(polygon, num.asarray(point)[num.newaxis, :])[0])

test/test_orthodrome.py

@@ -1,4 +1,5 @@
 from __future__ import division, print_function, absolute_import
+import os
 import unittest
 import math
 import random
@@ -22,7 +23,7 @@ r2d = 180./math.pi
 d2r = 1./r2d
 km = 1000.
 
-plot = False
+plot = int(os.environ.get('MPL_SHOW', 0))
 
 
 def random_lat(mi=-90., ma=90., rstate=None, size=None):
@@ -36,9 +37,18 @@ def random_lat(mi=-90., ma=90., rstate=None, size=None):
 def random_lon(mi=-180., ma=180., rstate=None, size=None):
     if rstate is None:
         rstate = num.random
-    mi_ = 0.5*(math.sin(mi * math.pi/180.)+1.)
-    ma_ = 0.5*(math.sin(ma * math.pi/180.)+1.)
-    return num.arcsin(rstate.uniform(mi_, ma_, size=size)*2.-1.)*180./math.pi
+
+    return rstate.uniform(mi, ma, size=size)
+
+
+def random_circle(npoints=100):
+    radius = num.random.uniform(0*km, 20000*km)
+    lat0, lon0 = random_lat(), random_lon()
+    phis = num.linspace(0., 2.*math.pi, npoints, endpoint=False)
+    ns, es = radius*num.sin(phis), radius*num.cos(phis)
+    lats, lons = orthodrome.ne_to_latlon(lat0, lon0, ns, es)
+    circle = num.vstack([lats, lons]).T
+    return lat0, lon0, radius, circle
 
 
 def light(color, factor=0.5):
@@ -389,71 +399,50 @@ class OrthodromeTestCase(unittest.TestCase):
         assert num.all(points2[:, 1] > -eps)
 
     def test_point_in_polygon(self):
-        if plot:
-            from pyrocko.plot import mpl_graph_color
 
+        if plot:
             import matplotlib.pyplot as plt
-            from matplotlib.patches import Polygon
-
-            axes = plt.gca()
-
-        nip = 100
-
-        for i in range(1):
-            np = 3
-            points = num.zeros((np, 2))
-            points[:, 0] = random_lat(size=3)
-            points[:, 1] = random_lon(size=3)
-
-            points_ip = num.zeros((nip*points.shape[0], 2))
-            for ip in range(points.shape[0]):
-                n, e = orthodrome.latlon_to_ne_numpy(
-                    points[ip % np, 0], points[ip % np, 1],
-                    points[(ip+1) % np, 0], points[(ip+1) % np, 1])
-
-                ns = num.arange(nip) * n / nip
-                es = num.arange(nip) * e / nip
-                lats, lons = orthodrome.ne_to_latlon(
-                    points[ip % np, 0], points[ip % np, 1], ns, es)
-
-                points_ip[ip*nip:(ip+1)*nip, 0] = lats
-                points_ip[ip*nip:(ip+1)*nip, 1] = lons
 
+        for i in range(100):
             if plot:
-                color = mpl_graph_color(i)
-                axes.add_patch(
-                    Polygon(
-                        num.fliplr(points_ip),
-                        facecolor=light(color),
-                        edgecolor=color,
-                        alpha=0.5))
+                plt.clf()
+                axes = plt.gca()
 
-            points_xyz = orthodrome.latlon_to_xyz(points_ip.T)
-            center_xyz = num.mean(points_xyz, axis=0)
+            lat0, lon0, radius, circle = random_circle(100)
+            if plot:
+                print(lat0, lon0, radius)
 
-            assert num.all(
-                orthodrome.distances3d(
-                    points_xyz, center_xyz[num.newaxis, :]) < 1.0)
+            lats = num.linspace(-90., 90., 100)
+            lons = num.linspace(-180., 180., 200)
 
-            lat, lon = orthodrome.xyz_to_latlon(center_xyz)
-            rot = orthodrome.rot_to_00(lat, lon)
+            points = num.empty((lons.size*lats.size, 2))
+            points[:, 0] = num.repeat(lats, lons.size)
+            points[:, 1] = num.tile(lons, lats.size)
 
-            points_rot_xyz = num.dot(rot, points_xyz.T).T
-            points_rot_pro = orthodrome.stereographic(points_rot_xyz)  # noqa
+            mask = orthodrome.contains_points(circle, points)
+            distances = orthodrome.distance_accurate50m_numpy(
+                lat0, lon0, points[:, 0], points[:, 1])
 
-            poly_xyz = orthodrome.latlon_to_xyz(points_ip)
-            poly_rot_xyz = num.dot(rot, poly_xyz.T).T
-            groups = orthodrome.spoly_cut([poly_rot_xyz], axis=0)
-            num.zeros(points.shape[0], dtype=num.int)
+            mask2 = distances < radius
+            mask3 = num.logical_and(
+                num.not_equal(mask2, mask),
+                num.abs(distances - radius) > radius / 100.)
 
             if plot:
-                for group in groups:
-                    for poly_rot_group_xyz in group:
+                axes.plot(
+                    circle[:, 1], circle[:, 0], 'o',
+                    ms=1, color='black')
+                axes.plot(
+                    points[mask, 1], points[mask, 0], 'o',
+                    ms=1, alpha=0.2, color='black')
+                axes.plot(
+                    points[mask3, 1], points[mask3, 0], 'o',
+                    ms=1, color='red')
 
-                        axes.set_xlim(-180., 180.)
-                        axes.set_ylim(-90., 90.)
+            if plot:
+                plt.show()
 
-                    plt.show()
+            assert not num.any(mask3)
 
     def test_point_in_region(self):
         testdata = [