source: anuga_core/source/anuga/utilities/test_polygon.py @ 6534

Last change on this file since 6534 was 6534, checked in by ole, 15 years ago

Added new function is_inside_triangle. This helps the fitting functions.

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