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

Last change on this file since 8879 was 8819, checked in by steve, 12 years ago

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