source: trunk/anuga_core/source/anuga/geometry/test_polygon.py @ 8032

Last change on this file since 8032 was 8032, checked in by habili, 14 years ago

Added tests for polygon_overlap and not_polygon_overlap

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