Merge branch 'master' of github.com:AIFDR/inasafe

safe/common/interpolation2d.py

@@ -46,7 +46,7 @@ def interpolate2d(x, y, Z, points, mode='linear', bounds_error=False):
             * 'linear' - bilinear interpolation using the four
                   nearest neighbours (default)
 
-        * bounds_error: Boolean flag. If True (default) a BoundsErorr exception
+        * bounds_error: Boolean flag. If True (default) a BoundsError exception
               will be raised when interpolated values are requested
               outside the domain of the input data. If False, nan
               is returned for those values
@@ -58,7 +58,8 @@ def interpolate2d(x, y, Z, points, mode='linear', bounds_error=False):
 
     Notes:
         Input coordinates x and y are assumed to be monotonically increasing,
-        but need not be equidistantly spaced.
+        but need not be equidistantly spaced. No such assumption regarding
+        ordering of points is made.
 
         Z is assumed to have dimension M x N, where M = len(x) and N = len(y).
         In other words it is assumed that the x values follow the first
@@ -241,24 +242,33 @@ def check_inputs(x, y, Z, points, mode, bounds_error):
     eta = points[:, 1]
 
     if bounds_error:
-        msg = ('Interpolation point %f was less than the smallest value in '
-               'domain %f and bounds_error was requested.' % (xi[0], x[0]))
-        if xi[0] < x[0]:
+        xi0 = min(xi)
+        xi1 = max(xi)
+        eta0 = min(eta)
+        eta1 = max(eta)
+
+        msg = ('Interpolation point xi=%f was less than the smallest '
+               'value in domain (x=%f) and bounds_error was requested.'
+               % (xi0, x[0]))
+        if xi0 < x[0]:
             raise BoundsError(msg)
 
-        msg = ('Interpolation point %f was greater than the largest value in '
-               'domain %f and bounds_error was requested.' % (xi[-1], x[-1]))
-        if xi[-1] > x[-1]:
+        msg = ('Interpolation point xi=%f was greater than the largest '
+               'value in domain (x=%f) and bounds_error was requested.'
+               % (xi1, x[-1]))
+        if xi1 > x[-1]:
             raise BoundsError(msg)
 
-        msg = ('Interpolation point %f was less than the smallest value in '
-               'domain %f and bounds_error was requested.' % (eta[0], y[0]))
-        if eta[0] < y[0]:
+        msg = ('Interpolation point eta=%f was less than the smallest '
+               'value in domain (y=%f) and bounds_error was requested.'
+               % (eta0, y[0]))
+        if eta0 < y[0]:
             raise BoundsError(msg)
 
-        msg = ('Interpolation point %f was greater than the largest value in '
-               'domain %f and bounds_error was requested.' % (eta[-1], y[-1]))
-        if eta[-1] > y[-1]:
+        msg = ('Interpolation point eta=%f was greater than the largest '
+               'value in domain (y=%f) and bounds_error was requested.'
+               % (eta1, y[-1]))
+        if eta1 > y[-1]:
             raise BoundsError(msg)
 
     return x, y, Z, xi, eta

safe/common/polygon.py

@@ -53,14 +53,10 @@ def separate_points_by_polygon(points, polygon,
         * use_numpy: Use the fast numpy implementation
 
     Returns:
-        * indices: array of same length as points with indices of points
-              falling inside the polygon listed from the beginning and indices
-              of points falling outside listed from the end.
-
-        * count: count of points falling inside the polygon
-
-        The indices of points inside are obtained as indices[:count]
-        The indices of points outside are obtained as indices[count:]
+        * indices_inside_polygon: array of indices of points
+              falling inside the polygon
+        * indices_outside_polygon: array of indices of points
+              falling outside the polygon
 
     Raises: A generic Exception is raised for unexpected input.
 
@@ -148,23 +144,21 @@ def separate_points_by_polygon(points, polygon,
     inside_box = -outside_box
     candidate_points = points[inside_box]
 
-    # FIXME (Ole): I would like to return just indices_inside, indices_outside
-    # instead of the legacy of one array with a break point
-    # in the underlying _separate_by_points functions too
     if use_numpy:
-        indices, count = _separate_points_by_polygon(candidate_points,
-                                                     polygon,
-                                                     closed=closed)
+        func = _separate_points_by_polygon
     else:
-        indices, count = _separate_points_by_polygon_python(candidate_points,
-                                                            polygon,
-                                                            closed=closed)
+        func = _separate_points_by_polygon_python
 
-    # Map local indices from candidate points to global indices of all points
-    indices_inside_polygon = numpy.where(inside_box)[0][indices[:count]]
+    local_indices_inside, local_indices_outside = func(candidate_points,
+                                                       polygon,
+                                                       closed=closed)
 
+    # Map local indices from candidate points to global indices of all points
     indices_outside_box = numpy.where(outside_box)[0]
-    indices_in_box_outside_poly = numpy.where(inside_box)[0][indices[count:]]
+    indices_inside_box = numpy.where(inside_box)[0]
+
+    indices_inside_polygon = indices_inside_box[local_indices_inside]
+    indices_in_box_outside_poly = indices_inside_box[local_indices_outside]
     indices_outside_polygon = numpy.concatenate((indices_outside_box,
                                                  indices_in_box_outside_poly))
 
@@ -203,8 +197,8 @@ def _separate_points_by_polygon(points, polygon,
     N = polygon.shape[0]
     M = points.shape[0]
     if M == 0:
-        # If no points return 0-vector
-        return numpy.arange(0), 0
+        # If no points return two 0-vectors
+        return numpy.arange(0), numpy.arange(0)
 
     x = points[:, 0]
     y = points[:, 1]
@@ -262,7 +256,7 @@ def _separate_points_by_polygon(points, polygon,
     # Indices of outside points
     indices[inside_index:] = numpy.where(1 - inside)[0]
 
-    return indices, inside_index
+    return indices[:inside_index], indices[inside_index:]
 
 
 def _separate_points_by_polygon_python(points, polygon,
@@ -312,7 +306,7 @@ def _separate_points_by_polygon_python(points, polygon,
     inside_index = 0  # Keep track of points inside
     outside_index = M - 1  # Keep track of points outside (starting from end)
 
-    # Begin main loop (for each point) - FIXME (write as vector ops)
+    # Begin main loop (for each point)
     for k in range(M):
         x = points[k, 0]
         y = points[k, 1]
@@ -362,7 +356,7 @@ def _separate_points_by_polygon_python(points, polygon,
     indices[inside_index:] = tmp[::-1]
 
     # Return reference result
-    return indices, inside_index
+    return indices[:inside_index], indices[inside_index:]
 
 
 def point_on_line(points, line, rtol=1.0e-5, atol=1.0e-8):
@@ -740,8 +734,6 @@ def clip_line_by_polygon(line, polygon,
     #    * Calculate its midpoint
     #    * Determine if it is inside or outside clipping polygon
 
-    # FIXME (Ole): Vectorise
-
     # Loop through line segments
     inside_line_segments = []
     outside_line_segments = []
@@ -1034,7 +1026,7 @@ def intersection(line0, line1, rtol=1.0e-12, atol=1.0e-12):
                                        [x2, y2], [x3, y3])
         else:
             # Lines are parallel but aren't collinear
-            return 4, None  # FIXME (Ole): Add distance here instead of None
+            return 4, None
     else:
         # Lines are not parallel, check if they intersect
         u0 = u0 / denom

safe/common/test_interpolate.py

@@ -321,6 +321,7 @@ def test_linear_interpolation_outside_domain(self):
 
         # Try a range of combinations of points outside domain
         # with error_bounds True
+        print
         for lox in [x[0], x[0] - 1]:
             for hix in [x[-1], x[-1] + 1]:
                 for loy in [y[0], y[0] - 1]:
@@ -332,15 +333,15 @@ def test_linear_interpolation_outside_domain(self):
                         points = combine_coordinates(xis, etas)
 
                         if lox < x[0] or hix > x[-1] or \
-                                loy < x[0] or hiy > y[-1]:
+                                loy < y[0] or hiy > y[-1]:
                             try:
                                 vals = interpolate2d(x, y, A, points,
                                                      mode='linear',
                                                      bounds_error=True)
-                            except BoundsError:
-                                pass
+                            except BoundsError, e:
+                                assert 'bounds_error was requested' in str(e)
                             else:
-                                msg = 'Should have raise bounds error'
+                                msg = 'Should have raised bounds error'
                                 raise Exception(msg)
 
         # Try a range of combinations of points outside domain with

safe/common/test_polygon.py

@@ -199,10 +199,7 @@ def test_inside_polygon_main2(self):
         assert is_inside_polygon((0.5, 0.5), polygon)
         assert is_inside_polygon((10.5, 10.5), polygon)
         assert not is_inside_polygon((0, 5.5), polygon)
-
-        # FIXME: Fails if point is 5.5, 5.5 - maybe ok, but
-        # need to understand it
-        #assert is_inside_polygon((5.5, 5.5), polygon)
+        assert is_inside_polygon((5.5, 5.5), polygon)
 
         # Polygon with a hole
         polygon = [[-1, -1], [2, -1], [2, 2], [-1, 2], [-1, -1],
@@ -250,6 +247,20 @@ def test_inside_polygon_vector_version(self):
         res = inside_polygon(points, polygon)
         assert numpy.allclose(res, [0, 1, 2])
 
+    def test_separate_points_by_polygon0(self):
+        """Points can be separated by polygon
+        """
+
+        # Now try the vector formulation returning indices
+        polygon = [[0, 0], [1, 0], [0.5, -1], [2, -1], [2, 1], [0, 1]]
+        points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+
+        inside, outside = separate_points_by_polygon(points, polygon)
+
+        assert len(inside) + len(outside) == len(points)
+        assert numpy.allclose(inside, [0, 1, 2])
+        assert numpy.allclose(outside, [3, 4])
+
     def test_outside_polygon(self):
         """Points are classified as either outside polygon or not
         """

safe/common/testing.py

@@ -4,6 +4,8 @@
 import numpy
 import os
 
+from numerics import axes2points
+
 # Find parent parent directory to path
 # NOTE: This must match Makefile target testdata
 # FIXME (Ole): Use environment variable for this.
@@ -24,12 +26,10 @@
 EXPDATA = os.path.join(DATADIR, 'exposure')  # Real exposure layers
 
 UNITDATA = os.path.abspath(
-                            #common/
-                            os.path.join(os.path.dirname(__file__),
-                            #safe/
-                            '..',
-                            'test',
-                            'data'))
+    os.path.join(os.path.dirname(__file__),
+                 '..',
+                 'test',
+                 'data'))
 
 # Known feature counts in test data
 FEATURE_COUNTS = {'test_buildings.shp': 144,
@@ -51,16 +51,8 @@ def combine_coordinates(x, y):
     """Make list of all combinations of points for x and y coordinates
     """
 
-    # FIXME (Ole): Write this using numpy for issue #91 and use that routine
-    # instead of the below. We need something like the Kronecker product
-    # numpy.kron
-    points = []
-    for px in x:
-        for py in y:
-            points.append((px, py))
-    points = numpy.array(points)
+    return axes2points(x, y)
 
-    return points
 
 # For polygon testing
 test_lines = [numpy.array([[122.231021, -8.626557],