source: trunk/anuga_core/anuga/geometry/tests/test_polygon.py @ 9680

Last change on this file since 9680 was 9566, checked in by steve, 10 years ago

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