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

Last change on this file since 8038 was 8038, checked in by habili, 13 years ago

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