# source:trunk/anuga_core/source/anuga/geometry/test_polygon.py@7866

Last change on this file since 7866 was 7866, checked in by hudson, 13 years ago

More swb tests passing. Cleaned up some pylint errors.

File size: 66.9 KB
RevLine
[5897]1#!/usr/bin/env python
2
[7866]3""" Test suite to test polygon functionality. """
4
[5897]5import unittest
[7276]6import numpy as num
[5897]7from anuga.utilities.numerical_tools import ensure_numeric
8from anuga.utilities.system_tools import get_pathname_from_package
9
[7866]10from polygon import _poly_xy, separate_points_by_polygon, \
11                    populate_polygon, polygon_area, is_inside_polygon, \
13                    plot_polygons, outside_polygon, is_outside_polygon, \
14                    intersection, is_complex, is_inside_triangle, \
15                    interpolate_polyline, inside_polygon, in_and_outside_polygon
16
[7858]17from polygon_function import Polygon_function
[5897]18from anuga.coordinate_transforms.geo_reference import Geo_reference
19from anuga.geospatial_data.geospatial_data import Geospatial_data
20
[7276]21
[5897]22def test_function(x, y):
23    return x+y
24
[7276]25
[5897]26class Test_Polygon(unittest.TestCase):
27    def setUp(self):
28        pass
29
30    def tearDown(self):
31        pass
32
33    def test_that_C_extension_compiles(self):
34        FN = 'polygon_ext.c'
35        try:
36            import polygon_ext
37        except:
38            from compile import compile
39
40            try:
41                compile(FN)
42            except:
[7276]43                raise Exception, 'Could not compile %s' %FN
[5897]44            else:
45                import polygon_ext
46
47    # Polygon stuff
48    def test_polygon_function_constants(self):
49        p1 = [[0,0], [10,0], [10,10], [0,10]]
50        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
51
[7276]52        f = Polygon_function([(p1, 1.0)])
53        z = f([5, 5, 27, 35], [5, 9, 8, -5])    # Two first inside p1
54        assert num.allclose(z, [1, 1, 0, 0])
[5897]55
[7276]56        f = Polygon_function([(p2, 2.0)])
57        z = f([5, 5, 27, 35], [5, 9, 8, -5])    # First and last inside p2
58        assert num.allclose(z, [2, 0, 0, 2])
[5897]59
[7276]60        # Combined
61        f = Polygon_function([(p1, 1.0), (p2, 2.0)])
62        z = f([5, 5, 27, 35], [5, 9, 8, -5])
63        assert num.allclose(z, [2, 1, 0, 2])
[5897]64
65    def test_polygon_function_csvfile(self):
66        from os import sep, getenv
67
68        # Get path where this test is run
69        path = get_pathname_from_package('anuga.utilities')
70
71        # Form absolute filename and read
[7276]72        filename = path + sep + 'mainland_only.csv'
[5897]74
[7276]75        f = Polygon_function([(p1, 10.0)])
76        z = f([430000, 480000], [490000, 7720000]) # 1st outside, 2nd inside
[5897]77
[7276]78        assert num.allclose(z, [0, 10])
79
[5897]80    def test_polygon_function_georef(self):
[7276]81        """Check that georeferencing works"""
[5897]82
83        from anuga.coordinate_transforms.geo_reference import Geo_reference
84
85        geo = Geo_reference(56, 200, 1000)
86
87        # Make points 'absolute'
88        p1 = [[200,1000], [210,1000], [210,1010], [200,1010]]
89        p2 = [[200,1000], [210,1010], [215,1005], [220, 1010], [225,1000],
90              [230,1010], [240,990]]
91
[7276]92        f = Polygon_function([(p1, 1.0)], geo_reference=geo)
93        z = f([5, 5, 27, 35], [5, 9, 8, -5]) # Two first inside p1
[5897]94
[7276]95        assert num.allclose(z, [1, 1, 0, 0])
[5897]96
[7276]97        f = Polygon_function([(p2, 2.0)], geo_reference=geo)
98        z = f([5, 5, 27, 35], [5, 9, 8, -5]) # First and last inside p2
99        assert num.allclose(z, [2, 0, 0, 2])
[5897]100
[7276]101        # Combined
102        f = Polygon_function([(p1, 1.0), (p2, 2.0)], geo_reference=geo)
103        z = f([5, 5, 27, 35], [5, 9, 8, -5])
104        assert num.allclose(z, [2, 1, 0, 2])
[5897]105
[7276]106        # Check that it would fail without geo
107        f = Polygon_function([(p1, 1.0), (p2, 2.0)])
108        z = f([5, 5, 27, 35], [5, 9, 8, -5])
109        assert not num.allclose(z, [2, 1, 0, 2])
[5897]110
[7276]111    def test_polygon_function_callable(self):
112        """Check that values passed into Polygon_function can be callable."""
[5897]113
114        p1 = [[0,0], [10,0], [10,10], [0,10]]
115        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
116
[7276]117        f = Polygon_function([(p1, test_function)])
118        z = f([5, 5, 27, 35], [5, 9, 8, -5]) # Two first inside p1
119        assert num.allclose(z, [10, 14, 0, 0])
[5897]120
[7276]121        # Combined
122        f = Polygon_function([(p1, test_function), (p2, 2.0)])
123        z = f([5, 5, 27, 35], [5, 9, 8, -5])
124        assert num.allclose(z, [2, 14, 0, 2])
[5897]125
[7276]126        # Combined w default
127        f = Polygon_function([(p1, test_function), (p2, 2.0)], default = 3.14)
128        z = f([5, 5, 27, 35], [5, 9, 8, -5])
129        assert num.allclose(z, [2, 14, 3.14, 2])
[5897]130
[7276]131        # Combined w default func
132        f = Polygon_function([(p1, test_function), (p2, 2.0)],
133                             default=test_function)
134        z = f([5, 5, 27, 35], [5, 9, 8, -5])
135        assert num.allclose(z, [2, 14, 35, 2])
[5897]136
137    def test_point_on_line(self):
[7276]138        # Endpoints first
139        assert point_on_line([0,0], [[0,0], [1,0]])
140        assert point_on_line([1,0], [[0,0], [1,0]])
[5897]141
[7276]142        # Endpoints first - non-simple
143        assert point_on_line([-0.1,0.0], [[-0.1,0.0], [1.5,0.6]])
144        assert point_on_line([1.5,0.6], [[-0.1,0.0], [1.5,0.6]])
[5897]145
[7276]146        # Then points on line
147        assert point_on_line([0.5,0], [[0,0], [1,0]])
148        assert point_on_line([0,0.5], [[0,1], [0,0]])
149        assert point_on_line([1,0.5], [[1,1], [1,0]])
150        assert point_on_line([0.5,0.5], [[0,0], [1,1]])
[5897]151
[7276]152        # Then points not on line
153        assert not point_on_line([0.5,0], [[0,1], [1,1]])
154        assert not point_on_line([0,0.5], [[0,0], [1,1]])
[5897]155
[7276]156        # From real example that failed
157        assert not point_on_line([40,50], [[40,20], [40,40]])
[5897]158
[7276]159        # From real example that failed
160        assert not point_on_line([40,19], [[40,20], [40,40]])
[5897]161
162        # Degenerate line
[7276]163        assert point_on_line([40,19], [[40,19], [40,19]])
164        assert not point_on_line([40,19], [[40,40], [40,40]])
[5897]165
166    def test_is_inside_polygon_main(self):
167        # Simplest case: Polygon is the unit square
168        polygon = [[0,0], [1,0], [1,1], [0,1]]
169
[7276]170        assert is_inside_polygon((0.5, 0.5), polygon)
171        assert not is_inside_polygon((0.5, 1.5), polygon)
172        assert not is_inside_polygon((0.5, -0.5), polygon)
173        assert not is_inside_polygon((-0.5, 0.5), polygon)
174        assert not is_inside_polygon((1.5, 0.5), polygon)
[5897]175
[7276]176        # Try point on borders
177        assert is_inside_polygon((1., 0.5), polygon, closed=True)
178        assert is_inside_polygon((0.5, 1.), polygon, closed=True)
179        assert is_inside_polygon((0., 0.5), polygon, closed=True)
180        assert is_inside_polygon((0.5, 0.), polygon, closed=True)
[5897]181
[7276]182        assert not is_inside_polygon((0.5, 1.), polygon, closed=False)
183        assert not is_inside_polygon((0., 0.5), polygon, closed=False)
184        assert not is_inside_polygon((0.5, 0.), polygon, closed=False)
185        assert not is_inside_polygon((1., 0.5), polygon, closed=False)
[5897]186
[7276]187    def test_inside_polygon_main(self):
[7841]188        """test_is_inside_polygon
[6534]189
190        Test fast version of of is_inside_polygon
191        """
192
193        # Simplest case: Polygon is the unit square
194        polygon = [[0,0], [1,0], [1,1], [0,1]]
195
[7841]196        assert is_inside_polygon( (0.5, 0.5), polygon )
197        assert not is_inside_polygon( (0.5, 1.5), polygon )
198        assert not is_inside_polygon( (0.5, -0.5), polygon )
199        assert not is_inside_polygon( (-0.5, 0.5), polygon )
200        assert not is_inside_polygon( (1.5, 0.5), polygon )
[6534]201
[7276]202        # Try point on borders
[7841]203        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
204        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
205        assert is_inside_polygon( (0., 0.5), polygon, closed=True)
206        assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
[6534]207
[7841]208        assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
209        assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
210        assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
211        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
[6534]212
213
[5897]214    def test_inside_polygon_main(self):
215        # Simplest case: Polygon is the unit square
[7276]216        polygon = [[0,0], [1,0], [1,1], [0,1]]
[5897]217
218        # From real example (that failed)
[7276]219        polygon = [[20,20], [40,20], [40,40], [20,40]]
220        points = [[40, 50]]
221        res = inside_polygon(points, polygon)
222        assert len(res) == 0
[5897]223
[7276]224        polygon = [[20,20], [40,20], [40,40], [20,40]]
225        points = [[25, 25], [30, 20], [40, 50], [90, 20], [40, 90]]
226        res = inside_polygon(points, polygon)
227        assert len(res) == 2
228        assert num.allclose(res, [0,1])
[5897]229
[7276]230        # More convoluted and non convex polygon
231        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
232        assert is_inside_polygon( (0.5, 0.5), polygon )
233        assert is_inside_polygon( (1, -0.5), polygon )
234        assert is_inside_polygon( (1.5, 0), polygon )
[5897]235
[7276]236        assert not is_inside_polygon( (0.5, 1.5), polygon )
237        assert not is_inside_polygon( (0.5, -0.5), polygon )
[5897]238
[7841]239        assert is_inside_polygon( (0.5, 0.5), polygon )
240        assert is_inside_polygon( (1, -0.5), polygon )
241        assert is_inside_polygon( (1.5, 0), polygon )
[5897]242
[7841]243        assert not is_inside_polygon( (0.5, 1.5), polygon )
244        assert not is_inside_polygon( (0.5, -0.5), polygon )
[5897]245
[7276]246        # Very convoluted polygon
[5897]247        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
[7276]248        assert is_inside_polygon((5, 5), polygon)
249        assert is_inside_polygon((17, 7), polygon)
250        assert is_inside_polygon((27, 2), polygon)
251        assert is_inside_polygon((35, -5), polygon)
252        assert not is_inside_polygon((15, 7), polygon)
253        assert not is_inside_polygon((24, 3), polygon)
254        assert not is_inside_polygon((25, -10), polygon)
[5897]255
[7276]256        # Another combination (that failed)
[5897]257        polygon = [[0,0], [10,0], [10,10], [0,10]]
[7276]258        assert is_inside_polygon((5, 5), polygon)
259        assert is_inside_polygon((7, 7), polygon)
260        assert not is_inside_polygon((-17, 7), polygon)
261        assert not is_inside_polygon((7, 17), polygon)
262        assert not is_inside_polygon((17, 7), polygon)
263        assert not is_inside_polygon((27, 8), polygon)
264        assert not is_inside_polygon((35, -5), polygon)
[5897]265
[7276]266        # Multiple polygons
267        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0], [10,10],
268                   [11,10], [11,11], [10,11], [10,10]]
269        assert is_inside_polygon((0.5, 0.5), polygon)
270        assert is_inside_polygon((10.5, 10.5), polygon)
[5897]271
[7276]272        # FIXME: Fails if point is 5.5, 5.5
273        assert not is_inside_polygon((0, 5.5), polygon)
[5897]274
[7276]275        # Polygon with a hole
[5897]276        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
[7276]277                   [0,0], [1,0], [1,1], [0,1], [0,0]]
[5897]278
[7276]279        assert is_inside_polygon((0, -0.5), polygon)
280        assert not is_inside_polygon((0.5, 0.5), polygon)
[5897]281
282    def test_duplicate_points_being_ok(self):
283        # Simplest case: Polygon is the unit square
284        polygon = [[0,0], [1,0], [1,0], [1,0], [1,1], [0,1], [0,0]]
285
[7276]286        assert is_inside_polygon((0.5, 0.5), polygon)
287        assert not is_inside_polygon((0.5, 1.5), polygon)
288        assert not is_inside_polygon((0.5, -0.5), polygon)
289        assert not is_inside_polygon((-0.5, 0.5), polygon)
290        assert not is_inside_polygon((1.5, 0.5), polygon)
[5897]291
[7276]292        # Try point on borders
293        assert is_inside_polygon((1., 0.5), polygon, closed=True)
294        assert is_inside_polygon((0.5, 1), polygon, closed=True)
295        assert is_inside_polygon((0., 0.5), polygon, closed=True)
296        assert is_inside_polygon((0.5, 0.), polygon, closed=True)
[5897]297
[7276]298        assert not is_inside_polygon((0.5, 1), polygon, closed=False)
299        assert not is_inside_polygon((0., 0.5), polygon, closed=False)
300        assert not is_inside_polygon((0.5, 0.), polygon, closed=False)
301        assert not is_inside_polygon((1., 0.5), polygon, closed=False)
[5897]302
303        # From real example
[7276]304        polygon = [[20,20], [40,20], [40,40], [20,40]]
305        points = [[40, 50]]
306        res = inside_polygon(points, polygon)
307        assert len(res) == 0
[5897]308
309    def test_inside_polygon_vector_version(self):
[7276]310        # Now try the vector formulation returning indices
[5897]311        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[7276]312        points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
313        res = inside_polygon( points, polygon, verbose=False )
[7308]314        assert num.allclose(res, [0,1,2])
[5897]315
316    def test_outside_polygon(self):
[7276]317        # unit square
318        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]319
320        # evaluate to False as the point 0.5, 0.5 is inside the unit square
[7276]321        assert not is_outside_polygon([0.5, 0.5], U)
322
[5897]323        # evaluate to True as the point 1.5, 0.5 is outside the unit square
[7276]324        assert is_outside_polygon([1.5, 0.5], U)
325
326        indices = outside_polygon([[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
327        assert num.allclose(indices, [1])
328
[5897]329        # One more test of vector formulation returning indices
330        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[7276]331        points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
332        res = outside_polygon(points, polygon)
333        assert num.allclose(res, [3, 4])
[5897]334
[7276]335        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
336        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
337                  [0.5, 1.5], [0.5, -0.5]]
338        res = outside_polygon(points, polygon)
[5897]339
[7276]340        assert num.allclose(res, [0, 4, 5])
[5897]341
[7276]342    def test_outside_polygon2(self):
343        # unit square
344        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]345
[7276]346        # evaluate to False as the point 0.5, 1.0 is inside the unit square
347        assert not outside_polygon([0.5, 1.0], U, closed = True)
[5897]348
349        # evaluate to True as the point 0.5, 1.0 is outside the unit square
[7276]350        assert is_outside_polygon([0.5, 1.0], U, closed = False)
[5897]351
352    def test_all_outside_polygon(self):
[7276]353        """Test case where all points are outside poly"""
[5897]354
[7276]355        # unit square
356        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]357
[7276]358        points = [[2,2], [1,3], [-1,1], [0,2]]      # All outside
[5897]359
360        indices, count = separate_points_by_polygon(points, U)
[7276]361        assert count == 0                           # None inside
362        assert num.allclose(indices, [3, 2, 1, 0])
[5897]363
364        indices = outside_polygon(points, U, closed = True)
[7276]365        assert num.allclose(indices, [0, 1, 2, 3])
[5897]366
367        indices = inside_polygon(points, U, closed = True)
[7276]368        assert num.allclose(indices, [])
[5897]369
370    def test_all_inside_polygon(self):
[7276]371        """Test case where all points are inside poly"""
[5897]372
[7276]373        # unit square
374        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]375
[7276]376        points = [[0.5,0.5], [0.2,0.3], [0,0.5]]    # All inside (or on edge)
[5897]377
378        indices, count = separate_points_by_polygon(points, U)
[7276]379        assert count == 3       # All inside
380        assert num.allclose(indices, [0, 1, 2])
[5897]381
[7276]382        indices = outside_polygon(points, U, closed=True)
[6158]383        assert num.allclose(indices, [])
[5897]384
[7276]385        indices = inside_polygon(points, U, closed=True)
386        assert num.allclose(indices, [0, 1, 2])
[5897]387
[7276]388
[5897]389    def test_separate_points_by_polygon(self):
[7276]390        # unit square
391        U = [[0,0], [1,0], [1,1], [0,1]]
[5897]392
[7276]393        indices, count = separate_points_by_polygon([[0.5, 0.5],
394                                                     [1, -0.5],
395                                                     [0.3, 0.2]],
396                                                    U)
397        assert num.allclose( indices, [0, 2, 1] )
[5897]398        assert count == 2
[7276]399
400        # One more test of vector formulation returning indices
[5897]401        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[7276]402        points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
403        res, count = separate_points_by_polygon(points, polygon)
[5897]404
[7276]405        assert num.allclose(res, [0, 1, 2, 4, 3])
406        assert count == 3
[5897]407
408        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[7866]409        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], \
410                        [0.5, -0.5]]
[7276]411        res, count = separate_points_by_polygon( points, polygon )
[5897]412
[7276]413        assert num.allclose( res, [1,2,3,5,4,0] )
414        assert count == 3
[5897]415
[7276]416
[6534]417    def test_is_inside_triangle(self):
418        # Simplest case:
[7276]419        triangle = [[0, 0], [1, 0], [0.5, 1]]
[6534]420
[7276]421        assert is_inside_triangle((0.5, 0.5), triangle)
422        assert is_inside_triangle((0.9, 0.1), triangle)
423        assert not is_inside_triangle((0.5, 1.5), triangle)
424        assert not is_inside_triangle((0.5, -0.5), triangle)
425        assert not is_inside_triangle((-0.5, 0.5), triangle)
426        assert not is_inside_triangle((1.5, 0.5), triangle)
[6534]427
[7276]428        # Try point on borders
429        assert is_inside_triangle((0.5, 0), triangle, closed=True)
430        assert is_inside_triangle((1, 0), triangle, closed=True)
[6534]431
[7276]432        assert not is_inside_triangle((0.5, 0), triangle, closed=False)
433        assert not is_inside_triangle((1, 0), triangle, closed=False)
[6534]434
435        # Try vertices
436        for P in triangle:
[7276]437            assert is_inside_triangle(P, triangle, closed=True)
438            assert not is_inside_triangle(P, triangle, closed=False)
439
[6534]440        # Slightly different
[7276]441        triangle = [[0, 0.1], [1, -0.2], [0.5, 1]]
442        assert is_inside_triangle((0.5, 0.5), triangle)
443        assert is_inside_triangle((0.4, 0.1), triangle)
444        assert not is_inside_triangle((1, 1), triangle)
445
[6534]446        # Try vertices
447        for P in triangle:
[7276]448            assert is_inside_triangle(P, triangle, closed=True)
449            assert not is_inside_triangle(P, triangle, closed=False)
[6534]450
[7276]451
[5897]452    def test_populate_polygon(self):
453        polygon = [[0,0], [1,0], [1,1], [0,1]]
454        points = populate_polygon(polygon, 5)
455
456        assert len(points) == 5
457        for point in points:
458            assert is_inside_polygon(point, polygon)
459
[7841]460            assert is_inside_polygon(point, polygon)
[5897]461
[7276]462
463        # Very convoluted polygon
[5897]464        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
465        points = populate_polygon(polygon, 5)
466        assert len(points) == 5
467        for point in points:
468            assert is_inside_polygon(point, polygon)
[7841]469            assert is_inside_polygon(point, polygon)
[5897]470
471
472    def test_populate_polygon_with_exclude(self):
473        polygon = [[0,0], [1,0], [1,1], [0,1]]
[7276]474        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]]     # SW quarter
475        points = populate_polygon(polygon, 5, exclude=[ex_poly])
[5897]476
477        assert len(points) == 5
478        for point in points:
479            assert is_inside_polygon(point, polygon)
[7276]480            assert not is_inside_polygon(point, ex_poly)
[5897]481
[7276]482        # overlap
[5897]483        polygon = [[0,0], [1,0], [1,1], [0,1]]
484        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
[7276]485        points = populate_polygon(polygon, 5, exclude=[ex_poly])
[5897]486
487        assert len(points) == 5
488        for point in points:
489            assert is_inside_polygon(point, polygon)
[7276]490            assert not is_inside_polygon(point, ex_poly)
491
492        # Multiple
[5897]493        polygon = [[0,0], [1,0], [1,1], [0,1]]
[7276]494        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]]    # SW quarter
495        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]]    # NE quarter
[5897]496
[7276]497        points = populate_polygon(polygon, 20, exclude=[ex_poly1, ex_poly2])
498
[5897]499        assert len(points) == 20
500        for point in points:
501            assert is_inside_polygon(point, polygon)
502            assert not is_inside_polygon(point, ex_poly1)
[7276]503            assert not is_inside_polygon(point, ex_poly2)
[5897]504
[7276]505        # Very convoluted polygon
[5897]506        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
507        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
[7276]508        points = populate_polygon(polygon, 20, exclude=[ex_poly])
509
[5897]510        assert len(points) == 20
511        for point in points:
512            assert is_inside_polygon(point, polygon)
[7276]513            assert not is_inside_polygon(point, ex_poly), '%s' % str(point)
[5897]514
515    def test_populate_polygon_with_exclude2(self):
[7276]516        min_outer = 0
[5897]517        max_outer = 1000
[7276]518        polygon_outer = [[min_outer,min_outer], [max_outer,min_outer],
519                         [max_outer,max_outer], [min_outer,max_outer]]
[5897]520
521        delta = 10
522        min_inner1 = min_outer + delta
523        max_inner1 = max_outer - delta
[7276]524        inner1_polygon = [[min_inner1,min_inner1], [max_inner1,min_inner1],
525                          [max_inner1,max_inner1], [min_inner1,max_inner1]]
[5897]526
[7276]527        density_inner2 = 1000
528        min_inner2 = min_outer + 2*delta
529        max_inner2 = max_outer - 2*delta
530        inner2_polygon = [[min_inner2,min_inner2], [max_inner2,min_inner2],
531                          [max_inner2,max_inner2], [min_inner2,max_inner2]]
532
533        points = populate_polygon(polygon_outer, 20,
534                                  exclude=[inner1_polygon, inner2_polygon])
535
[5897]536        assert len(points) == 20
537        for point in points:
538            assert is_inside_polygon(point, polygon_outer)
539            assert not is_inside_polygon(point, inner1_polygon)
[7276]540            assert not is_inside_polygon(point, inner2_polygon)
[5897]541
[7276]542        # Very convoluted polygon
[5897]543        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
544        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
[7276]545        points = populate_polygon(polygon, 20, exclude=[ex_poly])
546
[5897]547        assert len(points) == 20
548        for point in points:
549            assert is_inside_polygon(point, polygon)
[7276]550            assert not is_inside_polygon(point, ex_poly), '%s' % str(point)
[5897]551
552    def test_point_in_polygon(self):
553        polygon = [[0,0], [1,0], [1,1], [0,1]]
554        point = point_in_polygon(polygon)
555        assert is_inside_polygon(point, polygon)
556
[7276]557        # this may get into a vicious loop
558        # polygon = [[1e32,1e54], [1,0], [1,1], [0,1]]
559
[5897]560        polygon = [[1e15,1e7], [1,0], [1,1], [0,1]]
561        point = point_in_polygon(polygon)
562        assert is_inside_polygon(point, polygon)
563
564        polygon = [[0,0], [1,0], [1,1], [1e8,1e8]]
565        point = point_in_polygon(polygon)
566        assert is_inside_polygon(point, polygon)
567
568        polygon = [[1e32,1e54], [-1e32,1e54], [1e32,-1e54]]
569        point = point_in_polygon(polygon)
570        assert is_inside_polygon(point, polygon)
571
572        polygon = [[1e18,1e15], [1,0], [0,1]]
573        point = point_in_polygon(polygon)
574        assert is_inside_polygon(point, polygon)
575
576    def test_in_and_outside_polygon_main(self):
[7276]577        # Simplest case: Polygon is the unit square
578        polygon = [[0,0], [1,0], [1,1], [0,1]]
[5897]579
[7276]580        inside, outside = in_and_outside_polygon((0.5, 0.5), polygon)
581        assert inside[0] == 0
582        assert len(outside) == 0
[5897]583
[7276]584        inside, outside = in_and_outside_polygon((1., 0.5), polygon,
585                                                 closed=True)
586        assert inside[0] == 0
587        assert len(outside) == 0
[5897]588
[7276]589        inside, outside = in_and_outside_polygon((1., 0.5), polygon,
590                                                 closed=False)
591        assert len(inside) == 0
592        assert outside[0] == 0
[5897]593
[7276]594        points = [(1., 0.25),(1., 0.75)]
595        inside, outside = in_and_outside_polygon(points, polygon, closed=True)
[7308]596        assert num.alltrue(inside == [0,1])
[7276]597        assert len(outside) == 0
[5897]598
[7276]599        inside, outside = in_and_outside_polygon(points, polygon, closed=False)
600        assert len(inside) == 0
[7308]601        assert num.alltrue(outside == [0,1])
[5897]602
[7276]603        points = [(100., 0.25), (0.5, 0.5) ]
604        inside, outside = in_and_outside_polygon(points, polygon)
[7308]605        assert num.alltrue(inside == [1])
[7276]606        assert outside[0] == 0
[5897]607
[7276]608        points = [(100., 0.25),(0.5, 0.5), (39,20), (0.6,0.7),(56,43),(67,90)]
609        inside, outside = in_and_outside_polygon(points, polygon)
[7308]610        assert num.alltrue(inside == [1, 3])
611        assert num.alltrue(outside == [0, 2, 4, 5])
[7276]612
[5897]613    def test_intersection1(self):
614        line0 = [[-1,0], [1,0]]
615        line1 = [[0,-1], [0,1]]
616
617        status, value = intersection(line0, line1)
618        assert status == 1
[6158]619        assert num.allclose(value, [0.0, 0.0])
[5897]620
621    def test_intersection2(self):
622        line0 = [[0,0], [24,12]]
623        line1 = [[0,12], [24,0]]
624
625        status, value = intersection(line0, line1)
626        assert status == 1
[6158]627        assert num.allclose(value, [12.0, 6.0])
[5897]628
629        # Swap direction of one line
630        line1 = [[24,0], [0,12]]
631
632        status, value = intersection(line0, line1)
633        assert status == 1
[6158]634        assert num.allclose(value, [12.0, 6.0])
[5897]635
636        # Swap order of lines
637        status, value = intersection(line1, line0)
638        assert status == 1
[7276]639        assert num.allclose(value, [12.0, 6.0])
640
[5897]641    def test_intersection3(self):
642        line0 = [[0,0], [24,12]]
643        line1 = [[0,17], [24,0]]
644
645        status, value = intersection(line0, line1)
646        assert status == 1
[6158]647        assert num.allclose(value, [14.068965517, 7.0344827586])
[5897]648
649        # Swap direction of one line
650        line1 = [[24,0], [0,17]]
651
652        status, value = intersection(line0, line1)
653        assert status == 1
[7276]654        assert num.allclose(value, [14.068965517, 7.0344827586])
[5897]655
656        # Swap order of lines
657        status, value = intersection(line1, line0)
[7276]658        assert status == 1
659        assert num.allclose(value, [14.068965517, 7.0344827586])
[5897]660
661    def test_intersection_endpoints(self):
662        """test_intersection_endpoints(self):
663
664        Test that coinciding endpoints are picked up
665        """
[7276]666
[5897]667        line0 = [[0,0], [1,1]]
668        line1 = [[1,1], [2,1]]
669
670        status, value = intersection(line0, line1)
671        assert status == 1
[6158]672        assert num.allclose(value, [1.0, 1.0])
[5897]673
674        line0 = [[1,1], [2,0]]
675        line1 = [[1,1], [2,1]]
676
677        status, value = intersection(line0, line1)
678        assert status == 1
[7276]679        assert num.allclose(value, [1.0, 1.0])
[5897]680
[5942]681    # This function is a helper function for
682    # the test_intersection_bug_20081110_?() set of tests.
683    # This function tests all parallel line cases for 4 collinear points.
[5933]684    # This function should never be run directly by the unittest code.
685    def helper_test_parallel_intersection_code(self, P1, P2, P3, P4):
686        # lines in same direction, no overlap
687        # 0:         ---->----
688        # 1:                     --------->-----------
[7276]689        line0 = [P1, P2]
690        line1 = [P3, P4]
[5933]691        status, value = intersection(line0, line1)
692        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
693                               (str(status), str(value)))
694        self.failUnless(value is None, 'Expected value of None, got %s' %
695                                       str(value))
696
697        # lines in same direction, no overlap
698        # 0:         ----<----
699        # 1:                     ---------<-----------
[7276]700        line0 = [P2, P1]
701        line1 = [P4, P3]
[5933]702        status, value = intersection(line0, line1)
703        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
704                               (str(status), str(value)))
705        self.failUnless(value is None, 'Expected value of None, got %s' %
706                                       str(value))
707
708        # lines in opposite direction, no overlap
709        # 0:         ----<----
710        # 1:                     --------->-----------
[7276]711        line0 = [P2, P1]
712        line1 = [P3, P4]
[5933]713        status, value = intersection(line0, line1)
714        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
715                               (str(status), str(value)))
716        self.failUnless(value is None, 'Expected value of None, got %s' %
717                                       str(value))
718
719        # lines in opposite direction, no overlap
720        # 0:         ---->----
721        # 1:                     ---------<-----------
[7276]722        line0 = [P1, P2]
723        line1 = [P4, P3]
[5933]724        status, value = intersection(line0, line1)
725        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
726                               (str(status), str(value)))
727        self.failUnless(value is None, 'Expected value of None, got %s' %
728                                       str(value))
729
[7276]730        # ----------------------------------------------------------------------
[5942]731
[5933]732        # line0 fully within line1, same direction
733        # 0:         ---->----
734        # 1:    --------->-----------
735        # value should be line0:
736        #            ---->----
[7276]737        line0 = [P2, P3]
738        line1 = [P1, P4]
[5933]739        status, value = intersection(line0, line1)
740        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
741                               (str(status), str(value)))
[6158]742        self.failUnless(num.allclose(value, line0))
[5933]743
744        # line0 fully within line1, same direction
745        # 0:         ----<----
746        # 1:    ---------<-----------
747        # value should be line0:
748        #            ----<----
[7276]749        line0 = [P3, P2]
750        line1 = [P4, P1]
[5933]751        status, value = intersection(line0, line1)
752        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
753                               (str(status), str(value)))
[6158]754        self.failUnless(num.allclose(value, line0))
[5933]755
756        # line0 fully within line1, opposite direction
757        # 0:         ----<----
758        # 1:    --------->-----------
759        # value should be line0:
760        #            ----<----
[7276]761        line0 = [P3, P2]
762        line1 = [P1, P4]
[5933]763        status, value = intersection(line0, line1)
764        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
765                               (str(status), str(value)))
[6158]766        self.failUnless(num.allclose(value, line0))
[5933]767
768        # line0 fully within line1, opposite direction
769        # 0:         ---->----
770        # 1:    ---------<-----------
771        # value should be line0:
772        #            ---->----
[7276]773        line0 = [P2, P3]
774        line1 = [P4, P1]
[5933]775        status, value = intersection(line0, line1)
776        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
777                               (str(status), str(value)))
[6158]778        self.failUnless(num.allclose(value, line0))
[5933]779
[7276]780        # ----------------------------------------------------------------------
[5942]781
[5933]782        # line1 fully within line0, same direction
783        # 0:    --------->-----------
784        # 1:         ---->----
785        # value should be line1:
786        #            ---->----
[7276]787        line0 = [P1, P4]
788        line1 = [P2, P3]
[5933]789        status, value = intersection(line0, line1)
790        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
791                               (str(status), str(value)))
[6158]792        self.failUnless(num.allclose(value, line1))
[5933]793
794        # line1 fully within line0, same direction
795        # 0:    ---------<-----------
796        # 1:         ----<----
797        # value should be line1:
798        #            ----<----
[7276]799        line0 = [P4, P1]
800        line1 = [P3, P2]
[5933]801        status, value = intersection(line0, line1)
802        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
803                               (str(status), str(value)))
[6158]804        self.failUnless(num.allclose(value, line1))
[5933]805
806        # line1 fully within line0, opposite direction
807        # 0:    ---------<-----------
808        # 1:         ---->----
809        # value should be line1:
810        #            ---->----
[7276]811        line0 = [P4, P1]
812        line1 = [P2, P3]
[5933]813        status, value = intersection(line0, line1)
814        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
815                               (str(status), str(value)))
[6158]816        self.failUnless(num.allclose(value, line1))
[5933]817
818        # line1 fully within line0, opposite direction
819        # 0:    --------->-----------
820        # 1:         ----<----
821        # value should be line1:
822        #            ----<----
[7276]823        line0 = [P1, P4]
824        line1 = [P3, P2]
[5933]825        status, value = intersection(line0, line1)
826        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
827                               (str(status), str(value)))
[6158]828        self.failUnless(num.allclose(value, line1))
[5933]829
[7276]830        # ----------------------------------------------------------------------
[5942]831
[5933]832        # line in same direction, partial overlap
833        # 0:    ----->-----
834        # 1:       ------->--------
835        # value should be segment line1_start->line0_end:
836        #          --->----
[7276]837        line0 = [P1, P3]
838        line1 = [P2, P4]
[5933]839        status, value = intersection(line0, line1)
840        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
841                               (str(status), str(value)))
[6158]842        self.failUnless(num.allclose(value, [line1[0],line0[1]]))
[5933]843
844        # line in same direction, partial overlap
845        # 0:    -----<-----
846        # 1:       -------<--------
847        # value should be segment line0_start->line1_end:
848        #          ---<----
[7276]849        line0 = [P3, P1]
850        line1 = [P4, P2]
[5933]851        status, value = intersection(line0, line1)
852        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
853                               (str(status), str(value)))
[6158]854        self.failUnless(num.allclose(value, [line0[0],line1[1]]))
[5933]855
856        # line in opposite direction, partial overlap
857        # 0:    -----<-----
858        # 1:       ------->--------
859        # value should be segment line0_start->line1_start:
860        #          ---<----
[7276]861        line0 = [P3, P1]
862        line1 = [P2, P4]
[5933]863        status, value = intersection(line0, line1)
864        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
865                               (str(status), str(value)))
[6158]866        self.failUnless(num.allclose(value, [line0[0],line1[0]]))
[5933]867
868        # line in opposite direction, partial overlap
869        # 0:    ----->-----
870        # 1:       -------<--------
871        # value should be segment line1_end->line0_end:
872        #          --->----
[7276]873        line0 = [P1, P3]
874        line1 = [P4, P2]
[5933]875        status, value = intersection(line0, line1)
876        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
877                               (str(status), str(value)))
[6158]878        self.failUnless(num.allclose(value, [line1[1],line0[1]]))
[5933]879
[7276]880        # ----------------------------------------------------------------------
[5942]881
[5933]882        # line in same direction, partial overlap
883        # 0:       ------>------
884        # 1:    ------>------
885        # value should be segment line0_start->line1_end:
886        #          --->----
[7276]887        line0 = [P2, P4]
888        line1 = [P1, P3]
[5933]889        status, value = intersection(line0, line1)
890        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
891                               (str(status), str(value)))
[6158]892        self.failUnless(num.allclose(value, [line0[0],line1[1]]))
[5933]893
894        # line in same direction, partial overlap
895        # 0:       ------<------
896        # 1:    ------<------
897        # value should be segment line1_start->line0_end:
898        #          ----<-----
[7276]899        line0 = [P4, P2]
900        line1 = [P3, P1]
[5933]901        status, value = intersection(line0, line1)
902        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
903                               (str(status), str(value)))
[6158]904        self.failUnless(num.allclose(value, [line1[0],line0[1]]))
[5933]905
906        # line in opposite direction, partial overlap
907        # 0:       ------<------
908        # 1:    ----->------
909        # value should be segment line1_end->line0_end:
910        #          --->----
[7276]911        line0 = [P4, P2]
912        line1 = [P1, P3]
[5933]913        status, value = intersection(line0, line1)
914        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
915                               (str(status), str(value)))
[6158]916        self.failUnless(num.allclose(value, [line1[1],line0[1]]))
[5933]917
918        # line in opposite direction, partial overlap
919        # 0:       ------>------
920        # 1:    -----<------
921        # value should be segment line0_start->line1_start:
922        #          ---<----
[7276]923        line0 = [P2, P4]
924        line1 = [P3, P1]
[5933]925        status, value = intersection(line0, line1)
926        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
927                               (str(status), str(value)))
[6158]928        self.failUnless(num.allclose(value, [line0[0],line1[0]]))
[5933]929
[7276]930        # ----------------------------------------------------------------------
[5942]931
932        # line in same direction, same left point, line1 longer
933        # 0:    ----->------
934        # 1:    ------->--------
935        # value should be line0:
936        #       ----->------
[7276]937        line0 = [P1, P3]
938        line1 = [P1, P4]
[5942]939        status, value = intersection(line0, line1)
940        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
941                               (str(status), str(value)))
[6158]942        self.failUnless(num.allclose(value, line0))
[5942]943
944        # line in same direction, same left point, line1 longer
945        # 0:    -----<------
946        # 1:    -------<--------
947        # value should be line0:
948        #       -----<------
[7276]949        line0 = [P3, P1]
950        line1 = [P4, P1]
[5942]951        status, value = intersection(line0, line1)
952        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
953                               (str(status), str(value)))
[6158]954        self.failUnless(num.allclose(value, line0))
[5942]955
956        # line in opposite direction, same left point, line1 longer
957        # 0:    ----->------
958        # 1:    -------<--------
959        # value should be line0:
960        #       ----->------
[7276]961        line0 = [P1, P3]
962        line1 = [P4, P1]
[5942]963        status, value = intersection(line0, line1)
964        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
965                               (str(status), str(value)))
[6158]966        self.failUnless(num.allclose(value, line0))
[5942]967
968        # line in opposite direction, same start point, line1 longer
969        # 0:    -----<------
970        # 1:    ------->--------
971        # value should be line0:
972        #       -----<------
[7276]973        line0 = [P3, P1]
974        line1 = [P1, P4]
[5942]975        status, value = intersection(line0, line1)
976        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
977                               (str(status), str(value)))
[6158]978        self.failUnless(num.allclose(value, line0))
[5942]979
[7276]980        # ----------------------------------------------------------------------
[5942]981
982        # line in same direction, same left point, same right point
983        # 0:    ------->--------
984        # 1:    ------->--------
985        # value should be line0 or line1:
986        #       ------->--------
[7276]987        line0 = [P1, P3]
988        line1 = [P1, P3]
[5942]989        status, value = intersection(line0, line1)
990        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
991                               (str(status), str(value)))
[6158]992        self.failUnless(num.allclose(value, line0))
[5942]993
994        # line in same direction, same left point, same right point
995        # 0:    -------<--------
996        # 1:    -------<--------
997        # value should be line0 (or line1):
998        #       -------<--------
[7276]999        line0 = [P3, P1]
1000        line1 = [P3, P1]
[5942]1001        status, value = intersection(line0, line1)
1002        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1003                               (str(status), str(value)))
[6158]1004        self.failUnless(num.allclose(value, line0))
[5942]1005
1006        # line in opposite direction, same left point, same right point
1007        # 0:    ------->--------
1008        # 1:    -------<--------
1009        # value should be line0:
1010        #       ------->--------
[7276]1011        line0 = [P1, P3]
1012        line1 = [P3, P1]
[5942]1013        status, value = intersection(line0, line1)
1014        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1015                               (str(status), str(value)))
[6158]1016        self.failUnless(num.allclose(value, line0))
[5942]1017
1018        # line in opposite direction, same left point, same right point
1019        # 0:    -------<--------
1020        # 1:    ------->--------
1021        # value should be line0:
1022        #       -------<--------
[7276]1023        line0 = [P3, P1]
1024        line1 = [P1, P3]
[5942]1025        status, value = intersection(line0, line1)
1026        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1027                               (str(status), str(value)))
[6158]1028        self.failUnless(num.allclose(value, line0))
[5942]1029
[7276]1030        # ----------------------------------------------------------------------
[5942]1031
1032        # line in same direction, same right point, line1 longer
1033        # 0:        ----->------
1034        # 1:    ------->--------
1035        # value should be line0:
1036        #           ----->------
[7276]1037        line0 = [P2, P4]
1038        line1 = [P1, P4]
[5942]1039        status, value = intersection(line0, line1)
1040        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1041                               (str(status), str(value)))
[6158]1042        self.failUnless(num.allclose(value, line0))
[5942]1043
1044        # line in same direction, same right point, line1 longer
1045        # 0:        -----<------
1046        # 1:    -------<--------
1047        # value should be line0:
1048        #           -----<------
[7276]1049        line0 = [P4, P2]
1050        line1 = [P4, P1]
[5942]1051        status, value = intersection(line0, line1)
1052        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1053                               (str(status), str(value)))
[6158]1054        self.failUnless(num.allclose(value, line0))
[5942]1055
1056        # line in opposite direction, same right point, line1 longer
1057        # 0:        ----->------
1058        # 1:    -------<--------
1059        # value should be line0:
1060        #           ----->------
[7276]1061        line0 = [P2, P4]
1062        line1 = [P4, P1]
[5942]1063        status, value = intersection(line0, line1)
1064        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1065                               (str(status), str(value)))
[6158]1066        self.failUnless(num.allclose(value, line0))
[5942]1067
1068        # line in opposite direction, same right point, line1 longer
1069        # 0:        -----<------
1070        # 1:    ------->--------
1071        # value should be line0:
1072        #           -----<------
[7276]1073        line0 = [P4, P2]
1074        line1 = [P1, P4]
[5942]1075        status, value = intersection(line0, line1)
1076        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1077                               (str(status), str(value)))
[6158]1078        self.failUnless(num.allclose(value, line0))
[5942]1079
[7276]1080        # ----------------------------------------------------------------------
[5942]1081
1082        # line in same direction, same left point, line0 longer
1083        # 0:    ------->--------
1084        # 1:    ----->------
1085        # value should be line1:
1086        #       ----->------
[7276]1087        line0 = [P1, P4]
1088        line1 = [P1, P3]
[5942]1089        status, value = intersection(line0, line1)
1090        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1091                               (str(status), str(value)))
[6158]1092        self.failUnless(num.allclose(value, line1))
[5942]1093
1094        # line in same direction, same left point, line0 longer
1095        # 0:    -------<--------
1096        # 1:    -----<------
1097        # value should be line1:
1098        #       -----<------
[7276]1099        line0 = [P4, P1]
1100        line1 = [P3, P1]
[5942]1101        status, value = intersection(line0, line1)
1102        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1103                               (str(status), str(value)))
[6158]1104        self.failUnless(num.allclose(value, line1))
[5942]1105
1106        # line in opposite direction, same left point, line0 longer
1107        # 0:    ------->--------
1108        # 1:    -----<------
1109        # value should be line1:
1110        #       -----<------
[7276]1111        line0 = [P1, P4]
1112        line1 = [P3, P1]
[5942]1113        status, value = intersection(line0, line1)
1114        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1115                               (str(status), str(value)))
[6158]1116        self.failUnless(num.allclose(value, line1))
[5942]1117
1118        # line in opposite direction, same left point, line0 longer
1119        # 0:    -------<--------
1120        # 1:    ----->------
1121        # value should be line1:
1122        #       ----->------
[7276]1123        line0 = [P4, P1]
1124        line1 = [P1, P3]
[5942]1125        status, value = intersection(line0, line1)
1126        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1127                               (str(status), str(value)))
[6158]1128        self.failUnless(num.allclose(value, line1))
[5942]1129
[7276]1130        # ----------------------------------------------------------------------
[5942]1131
1132        # line in same direction, same right point, line0 longer
1133        # 0:    ------->--------
1134        # 1:        ----->------
1135        # value should be line1:
1136        #           ----->------
[7276]1137        line0 = [P1, P4]
1138        line1 = [P2, P4]
[5942]1139        status, value = intersection(line0, line1)
1140        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1141                               (str(status), str(value)))
[6158]1142        self.failUnless(num.allclose(value, line1))
[5942]1143
1144        # line in same direction, same right point, line0 longer
1145        # 0:    -------<--------
1146        # 1:        -----<------
1147        # value should be line1:
1148        #           -----<------
[7276]1149        line0 = [P4, P1]
1150        line1 = [P4, P2]
[5942]1151        status, value = intersection(line0, line1)
1152        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1153                               (str(status), str(value)))
[6158]1154        self.failUnless(num.allclose(value, line1))
[5942]1155
1156        # line in opposite direction, same right point, line0 longer
1157        # 0:    ------->--------
1158        # 1:        -----<------
1159        # value should be line1:
1160        #           -----<------
[7276]1161        line0 = [P1, P4]
1162        line1 = [P4, P2]
[5942]1163        status, value = intersection(line0, line1)
1164        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1165                               (str(status), str(value)))
[6158]1166        self.failUnless(num.allclose(value, line1))
[5942]1167
1168        # line in opposite direction, same right point, line0 longer
1169        # 0:    -------<--------
1170        # 1:        ----->------
1171        # value should be line1:
1172        #           ----->------
[7276]1173        line0 = [P4, P1]
1174        line1 = [P2, P4]
[5942]1175        status, value = intersection(line0, line1)
1176        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1177                               (str(status), str(value)))
[6158]1178        self.failUnless(num.allclose(value, line1))
[5942]1179
[7276]1180
[5942]1181    def test_intersection_bug_20081110_TR(self):
[5933]1182        """test_intersection_bug_20081110(self)
1183
[5942]1184        Test all cases in top-right quadrant
[5933]1185        """
1186
[5942]1187        # define 4 collinear points in top-right quadrant
[5933]1188        #    P1---P2---P3---P4
[5942]1189        P1 = [1.0, 1.0]
1190        P2 = [2.0, 2.0]
1191        P3 = [3.0, 3.0]
1192        P4 = [4.0, 4.0]
[7276]1193
[5942]1194        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1195        P1 = [1.0, 1.0+1.0e-9]
[7276]1196        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1197        P1 = [1.0, 1.0]
1198        P2 = [2.0, 2.0+1.0e-9]
[7276]1199        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1200        P2 = [2.0, 2.0]
1201        P3 = [3.0, 3.0+1.0e-9]
[7276]1202        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1203        P3 = [3.0, 3.0]
1204        P4 = [4.0, 4.0+1.0e-9]
[7276]1205        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1206
[5942]1207    def test_intersection_bug_20081110_TL(self):
1208        """test_intersection_bug_20081110(self)
[5933]1209
[5942]1210        Test all cases in top-left quadrant
1211        """
1212
1213        # define 4 collinear points in top-left quadrant
1214        #    P1---P2---P3---P4
1215        P1 = [-1.0, 1.0]
1216        P2 = [-2.0, 2.0]
1217        P3 = [-3.0, 3.0]
1218        P4 = [-4.0, 4.0]
[7276]1219
[5942]1220        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1221        P1 = [-1.0, 1.0+1.0e-9]
[7276]1222        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1223        P1 = [-1.0, 1.0]
1224        P2 = [-2.0, 2.0+1.0e-9]
[7276]1225        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1226        P2 = [-2.0, 2.0]
1227        P3 = [-3.0, 3.0+1.0e-9]
[7276]1228        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1229        P3 = [-3.0, 3.0]
1230        P4 = [-4.0, 4.0+1.0e-9]
[7276]1231        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1232
1233    def test_intersection_bug_20081110_BL(self):
[5933]1234        """test_intersection_bug_20081110(self)
1235
[5942]1236        Test all cases in bottom-left quadrant
[5933]1237        """
1238
[5942]1239        # define 4 collinear points in bottom-left quadrant
[5933]1240        #    P1---P2---P3---P4
[5942]1241        P1 = [-1.0, -1.0]
1242        P2 = [-2.0, -2.0]
1243        P3 = [-3.0, -3.0]
1244        P4 = [-4.0, -4.0]
[7276]1245
[5942]1246        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1247        P1 = [-1.0, -1.0+1.0e-9]
[7276]1248        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1249        P1 = [-1.0, -1.0]
1250        P2 = [-2.0, -2.0+1.0e-9]
[7276]1251        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1252        P2 = [-2.0, -2.0]
1253        P3 = [-3.0, -3.0+1.0e-9]
[7276]1254        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1255        P3 = [-3.0, -3.0]
1256        P4 = [-4.0, -4.0+1.0e-9]
[7276]1257        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1258
[5942]1259    def test_intersection_bug_20081110_BR(self):
[5933]1260        """test_intersection_bug_20081110(self)
1261
[5942]1262        Test all cases in bottom-right quadrant
[5933]1263        """
1264
[5942]1265        # define 4 collinear points in bottom-right quadrant
[5933]1266        #    P1---P2---P3---P4
[5942]1267        P1 = [1.0, -1.0]
1268        P2 = [2.0, -2.0]
1269        P3 = [3.0, -3.0]
1270        P4 = [4.0, -4.0]
[7276]1271
[5942]1272        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1273        P1 = [1.0, -1.0+1.0e-9]
[7276]1274        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1275        P1 = [1.0, -1.0]
1276        P2 = [2.0, -2.0+1.0e-9]
[7276]1277        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1278        P2 = [2.0, -2.0]
1279        P3 = [3.0, -3.0+1.0e-9]
[7276]1280        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5942]1281        P3 = [3.0, -3.0]
1282        P4 = [4.0, -4.0+1.0e-9]
[7276]1283        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1284
[5942]1285    def test_intersection_bug_20081110_TR_TL(self):
[5933]1286        """test_intersection_bug_20081110(self)
1287
[5942]1288        Test all cases in top-right & top-left quadrant
[5933]1289        """
1290
[5942]1291        # define 4 collinear points, 1 in TL, 3 in TR
[7276]1292        #    P1-+-P2---P3---P4
[5942]1293        P1 = [-3.0, 1.0]
1294        P2 = [ 1.0, 5.0]
1295        P3 = [ 2.0, 6.0]
1296        P4 = [ 3.0, 7.0]
1297        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1298
[5942]1299        # define 4 collinear points, 2 in TL, 2 in TR
[7276]1300        #    P1---P2-+-P3---P4
[5942]1301        P1 = [-3.0, 1.0]
1302        P2 = [-2.0, 2.0]
1303        P3 = [ 2.0, 6.0]
1304        P4 = [ 3.0, 7.0]
1305        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1306
[7276]1307        # define 4 collinear points, 3 in TL, 1 in TR
1308        #    P1---P2---P3-+-P4
[5942]1309        P1 = [-3.0, 1.0]
1310        P2 = [-2.0, 2.0]
1311        P3 = [-1.0, 3.0]
1312        P4 = [ 3.0, 7.0]
1313        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1314
[5942]1315    def test_intersection_bug_20081110_TR_BL(self):
[5933]1316        """test_intersection_bug_20081110(self)
1317
[5942]1318        Test all cases in top-right & bottom-left quadrant
[5933]1319        """
1320
[5942]1321        # define 4 collinear points, 1 in BL, 3 in TR
[7276]1322        #    P1-+-P2---P3---P4
[5942]1323        P1 = [-4.0, -3.0]
1324        P2 = [ 1.0,  2.0]
1325        P3 = [ 2.0,  3.0]
1326        P4 = [ 3.0,  4.0]
1327        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1328
[5942]1329        # define 4 collinear points, 2 in TL, 2 in TR
[7276]1330        #    P1---P2-+-P3---P4
[5942]1331        P1 = [-4.0, -3.0]
1332        P2 = [-3.0, -2.0]
1333        P3 = [ 2.0,  3.0]
1334        P4 = [ 3.0,  4.0]
1335        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1336
[5942]1337        # define 4 collinear points, 3 in TL, 1 in TR
[7276]1338        #    P1---P2---P3-+-P4
[5942]1339        P1 = [-4.0, -3.0]
1340        P2 = [-3.0, -2.0]
1341        P3 = [-2.0, -1.0]
1342        P4 = [ 3.0,  4.0]
1343        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5933]1344
[5942]1345    def test_intersection_bug_20081110_TR_BR(self):
1346        """test_intersection_bug_20081110(self)
[5939]1347
[5942]1348        Test all cases in top-right & bottom-right quadrant
[5939]1349        """
1350
[5942]1351        # define 4 collinear points, 1 in BR, 3 in TR
[7276]1352        #    P1-+-P2---P3---P4
[5942]1353        P1 = [ 1.0, -3.0]
1354        P2 = [ 5.0,  1.0]
1355        P3 = [ 6.0,  2.0]
1356        P4 = [ 7.0,  3.0]
1357        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5939]1358
[5942]1359        # define 4 collinear points, 2 in BR, 2 in TR
[7276]1360        #    P1---P2-+-P3---P4
[5942]1361        P1 = [ 1.0, -3.0]
1362        P2 = [ 2.0, -2.0]
1363        P3 = [ 6.0,  2.0]
1364        P4 = [ 7.0,  3.0]
1365        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5939]1366
[5942]1367        # define 4 collinear points, 3 in BR, 1 in TR
[7276]1368        #    P1---P2---P3-+-P4
[5942]1369        P1 = [ 1.0, -3.0]
1370        P2 = [ 2.0, -2.0]
1371        P3 = [ 3.0, -1.0]
1372        P4 = [ 7.0,  3.0]
1373        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
[5939]1374
[5940]1375
[5897]1376    def test_intersection_direction_invariance(self):
[7276]1377        """This runs through a number of examples and checks that
1378        direction of lines don't matter.
[5897]1379        """
[7276]1380
[5897]1381        line0 = [[0,0], [100,100]]
1382
1383        common_end_point = [20, 150]
[7276]1384
[5897]1385        for i in range(100):
1386            x = 20 + i * 1.0/100
1387
1388            line1 = [[x,0], common_end_point]
1389            status, p1 = intersection(line0, line1)
1390            assert status == 1
1391
1392            # Swap direction of line1
[7276]1393            line1 = [common_end_point, [x,0]]
[5897]1394            status, p2 = intersection(line0, line1)
[7276]1395            assert status == 1
[5897]1396
[7276]1397            msg = ('Orientation of line should not matter.\n'
1398                   'However, segment [%f,%f], [%f, %f]' %
1399                   (x, 0, common_end_point[0], common_end_point[1]))
1400            msg += ' gave %s, \nbut when reversed we got %s' % (p1, p2)
[6158]1401            assert num.allclose(p1, p2), msg
[5897]1402
1403            # Swap order of lines
1404            status, p3 = intersection(line1, line0)
[7276]1405            assert status == 1
[5897]1406            msg = 'Order of lines gave different results'
[6158]1407            assert num.allclose(p1, p3), msg
[5897]1408
1409    def test_no_intersection(self):
[7866]1410        """ Test 2 non-touching lines don't intersect. """
[5897]1411        line0 = [[-1,1], [1,1]]
1412        line1 = [[0,-1], [0,0]]
1413
1414        status, value = intersection(line0, line1)
1415        assert status == 0
1416        assert value is None
1417
1418    def test_intersection_parallel(self):
1419        line0 = [[-1,1], [1,1]]
1420        line1 = [[-1,0], [5,0]]
1421
1422        status, value = intersection(line0, line1)
[7276]1423        assert status == 4
[5897]1424        assert value is None
1425
1426        line0 = [[0,0], [10,100]]
1427        line1 = [[-10,5], [0,105]]
1428
1429        status, value = intersection(line0, line1)
[7276]1430        assert status == 4
1431        assert value is None
[5897]1432
1433    def test_intersection_coincide(self):
[7276]1434        """Test what happens when two lines partly coincide"""
[5897]1435
1436        # Overlap 1
1437        line0 = [[0,0], [5,0]]
1438        line1 = [[-3,0], [3,0]]
1439
1440        status, value = intersection(line0, line1)
1441        assert status == 2
[6158]1442        assert num.allclose(value, [[0,0], [3,0]])
[5897]1443
1444        # Overlap 2
1445        line0 = [[-10,0], [5,0]]
1446        line1 = [[-3,0], [10,0]]
1447
1448        status, value = intersection(line0, line1)
1449        assert status == 2
[7276]1450        assert num.allclose(value, [[-3, 0], [5,0]])
[5897]1451
1452        # Inclusion 1
1453        line0 = [[0,0], [5,0]]
1454        line1 = [[2,0], [3,0]]
1455
1456        status, value = intersection(line0, line1)
[7276]1457        assert status == 2
[6158]1458        assert num.allclose(value, line1)
[5897]1459
1460        # Inclusion 2
1461        line0 = [[1,0], [5,0]]
1462        line1 = [[-10,0], [15,0]]
1463
1464        status, value = intersection(line0, line1)
[7276]1465        assert status == 2
1466        assert num.allclose(value, line0)
[5897]1467
1468        # Exclusion (no intersection)
1469        line0 = [[-10,0], [1,0]]
1470        line1 = [[3,0], [15,0]]
1471
1472        status, value = intersection(line0, line1)
[7276]1473        assert status == 3
[5897]1474        assert value is None
1475
1476        # Try examples with some slope (y=2*x+5)
1477
1478        # Overlap
[7866]1479        line0 = [[0, 5], [7, 19]]
1480        line1 = [[1, 7], [10, 25]]
[5897]1481        status, value = intersection(line0, line1)
[7276]1482        assert status == 2
[6158]1483        assert num.allclose(value, [[1, 7], [7, 19]])
[5897]1484
1485        status, value = intersection(line1, line0)
1486        assert status == 2
[6158]1487        assert num.allclose(value, [[1, 7], [7, 19]])
[5897]1488
1489        # Swap direction
1490        line0 = [[7,19], [0,5]]
1491        line1 = [[1,7], [10,25]]
1492        status, value = intersection(line0, line1)
1493        assert status == 2
[6158]1494        assert num.allclose(value, [[7, 19], [1, 7]])
[5897]1495
1496        line0 = [[0,5], [7,19]]
1497        line1 = [[10,25], [1,7]]
1498        status, value = intersection(line0, line1)
1499        assert status == 2
[7276]1500        assert num.allclose(value, [[1, 7], [7, 19]])
[5897]1501
1502        # Inclusion
1503        line0 = [[1,7], [7,19]]
1504        line1 = [[0,5], [10,25]]
1505        status, value = intersection(line0, line1)
[7276]1506        assert status == 2
1507        assert num.allclose(value, [[1,7], [7, 19]])
[5897]1508
1509        line0 = [[0,5], [10,25]]
1510        line1 = [[1,7], [7,19]]
1511        status, value = intersection(line0, line1)
[7276]1512        assert status == 2
[6158]1513        assert num.allclose(value, [[1,7], [7, 19]])
[5897]1514
1515        line0 = [[0,5], [10,25]]
1516        line1 = [[7,19], [1,7]]
1517        status, value = intersection(line0, line1)
[7276]1518        assert status == 2
1519        assert num.allclose(value, [[7, 19], [1, 7]])
[5897]1520
1521
[7866]1522    def test_inside_polygon_main(self):
[7276]1523        # From real example (that failed)
1524        polygon = [[20,20], [40,20], [40,40], [20,40]]
1525        points = [[40, 50]]
1526        res = inside_polygon(points, polygon)
1527        assert len(res) == 0
[5897]1528
[7276]1529        polygon = [[20,20], [40,20], [40,40], [20,40]]
1530        points = [[25, 25], [30, 20], [40, 50], [90, 20], [40, 90]]
1531        res = inside_polygon(points, polygon)
1532        assert len(res) == 2
1533        assert num.allclose(res, [0,1])
[5897]1534
1535    def test_polygon_area(self):
[7866]1536        """ Test getting the area of a polygon. """
[7276]1537        # Simplest case: Polygon is the unit square
[5897]1538        polygon = [[0,0], [1,0], [1,1], [0,1]]
[7276]1539        assert polygon_area(polygon) == 1
[5897]1540
[7276]1541        # Simple case: Polygon is a rectangle
[5897]1542        polygon = [[0,0], [1,0], [1,4], [0,4]]
[7276]1543        assert polygon_area(polygon) == 4
[5897]1544
[7276]1545        # Simple case: Polygon is a unit triangle
[5897]1546        polygon = [[0,0], [1,0], [0,1]]
[7276]1547        assert polygon_area(polygon) == 0.5
[5897]1548
[7276]1549        # Simple case: Polygon is a diamond
[5897]1550        polygon = [[0,0], [1,1], [2,0], [1, -1]]
[7276]1551        assert polygon_area(polygon) == 2.0
1552
[6000]1553        # Complex case where numerical errors might occur
1554        polygon = [[314037.58727982, 6224952.2960092],
[7276]1555                   [314038.58727982, 6224952.2960092],
[6000]1556                   [314038.58727982, 6224953.2960092],
1557                   [314037.58727982, 6224953.2960092]]
[7276]1558        assert polygon_area(polygon) == 1.0
[5897]1559
1560    def test_poly_xy(self):
[7276]1561        # Simplest case: Polygon is the unit square
[5897]1562        polygon = [[0,0], [1,0], [1,1], [0,1]]
[7858]1563        x, y = _poly_xy(polygon)
[7276]1564        assert len(x) == len(polygon)+1
1565        assert len(y) == len(polygon)+1
1566        assert x[0] == 0
1567        assert x[1] == 1
1568        assert x[2] == 1
1569        assert x[3] == 0
1570        assert y[0] == 0
1571        assert y[1] == 0
1572        assert y[2] == 1
1573        assert y[3] == 1
[5897]1574
[7276]1575    # Arbitrary polygon
[5897]1576        polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]]
[7858]1577        x, y = _poly_xy(polygon)
[7276]1578        assert len(x) == len(polygon)+1
1579        assert len(y) == len(polygon)+1
1580        assert x[0] == 1
1581        assert x[1] == 1
1582        assert x[2] == 100
1583        assert x[3] == 1
1584        assert x[4] == 3
1585        assert y[0] == 5
1586        assert y[1] == 1
1587        assert y[2] == 10
1588        assert y[3] == 10
1589        assert y[4] == 6
[5897]1590
[7841]1591
1592    def test_plot_polygons(self):
[5897]1593        import os
[7276]1594
1595        # Simplest case: Polygon is the unit square
[5897]1596        polygon1 = [[0,0], [1,0], [1,1], [0,1]]
1597        polygon2 = [[1,1], [2,1], [3,2], [2,2]]
[7858]1598        plot_polygons([polygon1, polygon2], figname='test1')
[5897]1599
[7276]1600        # Another case
[5897]1601        polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
[7866]1602        plot_polygons([polygon2, polygon3], figname='test2')
[5897]1603
[7858]1604        for file in ['test1.png', 'test2.png']:
1605            assert os.access(file, os.R_OK)
1606            os.remove(file)
[5897]1607
[7858]1608
[5897]1609    def test_inside_polygon_geospatial(self):
[7866]1610        """ Test geospatial coords inside poly. """
[5897]1611        #Simplest case: Polygon is the unit square
[7866]1612        polygon_absolute = [[0, 0], [1, 0], [1, 1], [0, 1]]
[7276]1613        poly_geo_ref = Geo_reference(57, 100, 100)
[5897]1614        polygon = poly_geo_ref.change_points_geo_ref(polygon_absolute)
[7276]1615        poly_spatial = Geospatial_data(polygon, geo_reference=poly_geo_ref)
1616
[5897]1617        points_absolute = (0.5, 0.5)
[7276]1618        points_geo_ref = Geo_reference(57, 78, -56)
[5897]1619        points = points_geo_ref.change_points_geo_ref(points_absolute)
[7276]1620        points_spatial = Geospatial_data(points, geo_reference=points_geo_ref)
1621
[5897]1622        assert is_inside_polygon(points_absolute, polygon_absolute)
1623        assert is_inside_polygon(ensure_numeric(points_absolute),
1624                                 ensure_numeric(polygon_absolute))
[7276]1625        assert is_inside_polygon(points_absolute, poly_spatial)
1626        assert is_inside_polygon(points_spatial, poly_spatial)
1627        assert is_inside_polygon(points_spatial, polygon_absolute)
[5897]1628
[7276]1629        assert is_inside_polygon(points_absolute, polygon_absolute)
[5897]1630
[7866]1631    def test_decimate_polygon(self):
1632        from polygon import decimate_polygon
[5897]1633        polygon = [[0,0], [10,10], [15,5], [20, 10],
1634                   [25,0], [30,10], [40,-10], [35, -5]]
1635
1636        dpoly = decimate_polygon(polygon, factor=2)
1637
1638        assert len(dpoly)*2==len(polygon)
1639
[7276]1640
[6189]1641    def test_interpolate_polyline(self):
1642        """test_interpolate_polyline(self):
[7276]1643
1644        This test is added under the assumption that the function
1645        interpolate_polyline implemented by John Jakeman works.
1646        It has been exercised somewhat by tests of sts boundary,
1647        but never before separately.
[6189]1648        """
[7276]1649
[6189]1650        f = num.array([58.06150614, 58.06150614, 58.06150614])
1651        vertex_coordinates = num.array([[0., 0., ],
1652                                        [4.04092634, 1106.11074699],
1653                                        [8.08836552, 2212.16910609]])
1654        gauge_neighbour_id = [1, 2, -1]
1655        point_coordinates = num.array([[2.21870766e+03, 1.09802864e+03],
1656                                       [1.62739645e+03, 2.20626983e+03],
1657                                       [5.20084967e+02, 2.21030386e+03],
1658                                       [6.06464546e+00, 1.65913993e+03],
1659                                       [1.61934862e+03, -5.88143836e+00],
1660                                       [5.11996623e+02, -1.85956061e+00],
1661                                       [2.02046270e+00, 5.53055373e+02]])
[7276]1662
[6189]1663        z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
[7276]1664
1665        z = interpolate_polyline(f, vertex_coordinates,
1666                                 gauge_neighbour_id, point_coordinates)
[6189]1667        assert num.allclose(z, z_ref)
[7276]1668
[6189]1669        # Another f
1670        f = num.array([58.06150614, 158.06150614, 258.06150614])
[7276]1671        z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]
1672        z = interpolate_polyline(f, vertex_coordinates,
1673                                 gauge_neighbour_id, point_coordinates)
1674        assert num.allclose(z, z_ref)
[6189]1675
1676        # Other and simpler numbers
[7276]1677        f = num.array([1, 5, 13])
[6189]1678        vertex_coordinates = num.array([[0., 0., ],
1679                                        [4., 4.],
[7276]1680                                        [8., 8.]])
[6189]1681        point_coordinates = num.array([[0.1, 0.1],
1682                                       [3.5, 3.5],
1683                                       [4.0, 4.0],
1684                                       [5.2, 5.2],
1685                                       [7.0, 7.0],
1686                                       [8.3, 8.3]])
1687        gauge_neighbour_id = [1, 2, -1]
[7276]1688
1689        z = interpolate_polyline(f, vertex_coordinates,
1690                                 gauge_neighbour_id, point_coordinates)
[6189]1691        z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
[7276]1692        assert num.allclose(z, z_ref)
1693
[6189]1694        # Test exception thrown for one point
[7276]1695        f = num.array([5])
1696        vertex_coordinates = num.array([[4., 4.]])
1697        try:
1698            z = interpolate_polyline(f, vertex_coordinates,
1699                                     gauge_neighbour_id, point_coordinates)
[6189]1700        except Exception:
1701            pass
1702        else:
1703            raise Exception, 'One point should have raised exception'
1704
1705        # More polyline nodes
[7276]1706        data = num.array([1, 5, 13, 12, 6, 29])
[6189]1707        polyline_nodes = num.array([[0., 0.],
1708                                    [4., 4.],
1709                                    [8., 8.],
1710                                    [10., 10.],
1711                                    [10., 5.],
1712                                    [10., 0.]])
1713        point_coordinates = num.array([[0.1, 0.1],
1714                                       [3.5, 3.5],
1715                                       [4.0, 4.0],
1716                                       [5.2, 5.2],
1717                                       [7.0, 7.0],
1718                                       [8.3, 8.3],
1719                                       [10., 10.],
1720                                       [10., 9.],
1721                                       [10., 7.1],
1722                                       [10., 4.3],
1723                                       [10., 1.0]])
1724        gauge_neighbour_id = [1, 2, 3, 4, 5, -1]
[7276]1725        z = interpolate_polyline(data, polyline_nodes,
1726                                 gauge_neighbour_id, point_coordinates)
[6189]1727        z_ref = [1.1, 4.5, 5., 7.4, 11., 12.85, 12., 10.8, 8.52, 9.22, 24.4]
[7276]1728        assert num.allclose(z, z_ref)
[6535]1729
[7276]1730
1731    def test_is_inside_triangle_more(self):
[7866]1732        """ Test if points inside triangles are detected correctly. """
[6535]1733        res = is_inside_triangle([0.5, 0.5], [[ 0.5,  0. ],
1734                                              [ 0.5,  0.5],
1735                                              [ 0.,   0. ]])
1736        assert res is True
1737
1738        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1739                                 [[ 0.5,  0. ], [ 0.5, 0.5], [0., 0.]])
1740        assert res is False
[7276]1741
[6535]1742        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1743                                 [[1., 0.], [1., 0.5], [0.5, 0.]])
[7276]1744        assert res is False
1745
1746
[6535]1747        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1748                                 [[0.5, 0.5], [0.5, 0.], [1., 0.5]])
[7276]1749        assert res is True
[6535]1750
[7276]1751
[6535]1752        res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
1753                                 [[0.5, 0.], [0.5, 0.5], [0., 0.]])
1754        assert res is False
1755
[7276]1756
[6535]1757        res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
1758                                 [[0., 0.5], [0., 0.], [0.5, 0.5]])
[7276]1759        assert res is True
1760
[6535]1761        res = is_inside_triangle([0.69999999999999996, 0.69999999999999996],
1762                                 [[0.5, 0.], [0.5, 0.5], [0., 0.]])
[7276]1763        assert res is False
[6535]1764
1765        res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1766                                 [[0.25, 0.5], [0.25, 0.25], [0.5, 0.5]])
1767        assert res is False
[6544]1768
[7276]1769
1770        res = is_inside_triangle([10, 3],
[6544]1771                                 [[0.1, 0.],
1772                                  [0.1, 0.08333333],
1773                                  [0., 0.]])
1774        assert res is False
[7276]1775
[7687]1776
[7686]1777    def test_is_polygon_complex(self):
[7866]1778        """ Test a concave and a complex poly with is_complex, to make
1779            sure it can detect self-intersection.
[7690]1780        """
[7686]1781        concave_poly = [[0, 0], [10, 0], [5, 5], [10, 10], [0, 10]]
[7866]1782        complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], \
1783                            [0, 10]]
[7687]1784
[7690]1785        assert not is_complex(concave_poly)
1786        assert is_complex(complex_poly)
[7687]1787
1788    def test_is_polygon_complex2(self):
[7866]1789        """ Test a concave and a complex poly with is_complex, to make sure it
1790            can detect self-intersection. This test uses more complicated
1791            polygons.
[7690]1792        """
[7866]1793        concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7, 6], [10, 10], \
1794                        [1, 5], [0, 10]]
1795        complex_poly = [[0, 0], [12, 12], [10, 0], [5, 5], [3,18], [4, 15], \
1796                        [5, 7], [10, 10], [0, 10], [16, 12]]
[7687]1797
[7690]1798        assert not is_complex(concave_poly)
1799        assert is_complex(complex_poly)
[7686]1800
[7276]1801################################################################################
1802
[5897]1803if __name__ == "__main__":
[7276]1804    suite = unittest.makeSuite(Test_Polygon,'test')
[5897]1805    runner = unittest.TextTestRunner()
1806    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.