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

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

Work done during Water Down Under 2008.
Line Mesh intersections towards ticket:175 (flow through a cross section).
File size: 26.8 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
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 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 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 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 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 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 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 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 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 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 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 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 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
168
169
170    def test_is_inside_polygon_main(self):
171
172
173        # Simplest case: Polygon is the unit square
174        polygon = [[0,0], [1,0], [1,1], [0,1]]
175
176        assert is_inside_polygon( (0.5, 0.5), polygon )
177        assert not is_inside_polygon( (0.5, 1.5), polygon )
178        assert not is_inside_polygon( (0.5, -0.5), polygon )
179        assert not is_inside_polygon( (-0.5, 0.5), polygon )
180        assert not is_inside_polygon( (1.5, 0.5), polygon )
181
182        # Try point on borders
183        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
184        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
185        assert is_inside_polygon( (0., 0.5), polygon, closed=True)
186        assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
187
188        assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
189        assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
190        assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
191        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
192
193
194    def test_inside_polygon_main(self):
195
196        # Simplest case: Polygon is the unit square
197        polygon = [[0,0], [1,0], [1,1], [0,1]]       
198
199        # From real example (that failed)
200        polygon = [[20,20], [40,20], [40,40], [20,40]]
201        points = [ [40, 50] ]
202        res = inside_polygon(points, polygon)
203        assert len(res) == 0
204
205        polygon = [[20,20], [40,20], [40,40], [20,40]]
206        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
207        res = inside_polygon(points, polygon)
208        assert len(res) == 2
209        assert allclose(res, [0,1])
210
211
212
213        # More convoluted and non convex polygon
214        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
215        assert is_inside_polygon( (0.5, 0.5), polygon )
216        assert is_inside_polygon( (1, -0.5), polygon )
217        assert is_inside_polygon( (1.5, 0), polygon )
218
219        assert not is_inside_polygon( (0.5, 1.5), polygon )
220        assert not is_inside_polygon( (0.5, -0.5), polygon )
221
222
223        # Very convoluted polygon
224        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
225        assert is_inside_polygon( (5, 5), polygon )
226        assert is_inside_polygon( (17, 7), polygon )
227        assert is_inside_polygon( (27, 2), polygon )
228        assert is_inside_polygon( (35, -5), polygon )
229        assert not is_inside_polygon( (15, 7), polygon )
230        assert not is_inside_polygon( (24, 3), polygon )
231        assert not is_inside_polygon( (25, -10), polygon )
232
233
234
235        # Another combination (that failed)
236        polygon = [[0,0], [10,0], [10,10], [0,10]]
237        assert is_inside_polygon( (5, 5), polygon )
238        assert is_inside_polygon( (7, 7), polygon )
239        assert not is_inside_polygon( (-17, 7), polygon )
240        assert not is_inside_polygon( (7, 17), polygon )
241        assert not is_inside_polygon( (17, 7), polygon )
242        assert not is_inside_polygon( (27, 8), polygon )
243        assert not is_inside_polygon( (35, -5), polygon )
244
245
246
247
248        # Multiple polygons
249
250        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
251                   [10,10], [11,10], [11,11], [10,11], [10,10]]
252        assert is_inside_polygon( (0.5, 0.5), polygon )
253        assert is_inside_polygon( (10.5, 10.5), polygon )
254
255        #FIXME: Fails if point is 5.5, 5.5
256        assert not is_inside_polygon( (0, 5.5), polygon )
257
258        # Polygon with a hole
259        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
260                   [0,0], [1,0], [1,1], [0,1], [0,0]]
261
262        assert is_inside_polygon( (0, -0.5), polygon )
263        assert not is_inside_polygon( (0.5, 0.5), polygon )
264
265
266
267    def test_duplicate_points_being_ok(self):
268
269
270        # Simplest case: Polygon is the unit square
271        polygon = [[0,0], [1,0], [1,0], [1,0], [1,1], [0,1], [0,0]]
272
273        assert is_inside_polygon( (0.5, 0.5), polygon )
274        assert not is_inside_polygon( (0.5, 1.5), polygon )
275        assert not is_inside_polygon( (0.5, -0.5), polygon )
276        assert not is_inside_polygon( (-0.5, 0.5), polygon )
277        assert not is_inside_polygon( (1.5, 0.5), polygon )
278
279        # Try point on borders
280        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
281        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
282        assert is_inside_polygon( (0., 0.5), polygon, closed=True)
283        assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
284
285        assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
286        assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
287        assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
288        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
289
290        # From real example
291        polygon = [[20,20], [40,20], [40,40], [20,40]]
292        points = [ [40, 50] ]
293        res = inside_polygon(points, polygon)
294        assert len(res) == 0
295
296       
297
298    def test_inside_polygon_vector_version(self):
299        # Now try the vector formulation returning indices
300        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
301        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
302        res = inside_polygon( points, polygon, verbose=False )
303
304        assert allclose( res, [0,1,2] )
305
306    def test_outside_polygon(self):
307        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
308
309        assert not is_outside_polygon( [0.5, 0.5], U )
310        # evaluate to False as the point 0.5, 0.5 is inside the unit square
311       
312        assert is_outside_polygon( [1.5, 0.5], U )
313        # evaluate to True as the point 1.5, 0.5 is outside the unit square
314       
315        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
316        assert allclose( indices, [1] )
317       
318        # One more test of vector formulation returning indices
319        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
320        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
321        res = outside_polygon( points, polygon )
322
323        assert allclose( res, [3, 4] )
324
325
326
327        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
328        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
329        res = outside_polygon( points, polygon )
330
331        assert allclose( res, [0, 4, 5] )       
332     
333    def test_outside_polygon2(self):
334        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
335   
336        assert not outside_polygon( [0.5, 1.0], U, closed = True )
337        # evaluate to False as the point 0.5, 1.0 is inside the unit square
338       
339        assert is_outside_polygon( [0.5, 1.0], U, closed = False )
340        # evaluate to True as the point 0.5, 1.0 is outside the unit square
341
342    def test_all_outside_polygon(self):
343        """Test case where all points are outside poly
344        """
345       
346        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
347
348        points = [[2,2], [1,3], [-1,1], [0,2]] #All outside
349
350
351        indices, count = separate_points_by_polygon(points, U)
352        #print indices, count
353        assert count == 0 #None inside
354        assert allclose(indices, [3,2,1,0])
355
356        indices = outside_polygon(points, U, closed = True)
357        assert allclose(indices, [0,1,2,3])
358
359        indices = inside_polygon(points, U, closed = True)
360        assert allclose(indices, [])               
361
362
363    def test_all_inside_polygon(self):
364        """Test case where all points are inside poly
365        """
366       
367        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
368
369        points = [[0.5,0.5], [0.2,0.3], [0,0.5]] #All inside (or on edge)
370
371
372        indices, count = separate_points_by_polygon(points, U)
373        assert count == 3 #All inside
374        assert allclose(indices, [0,1,2])
375
376        indices = outside_polygon(points, U, closed = True)
377        assert allclose(indices, [])
378
379        indices = inside_polygon(points, U, closed = True)
380        assert allclose(indices, [0,1,2])
381       
382
383    def test_separate_points_by_polygon(self):
384        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
385
386        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
387        assert allclose( indices, [0,2,1] )
388        assert count == 2
389       
390        #One more test of vector formulation returning indices
391        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
392        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
393        res, count = separate_points_by_polygon( points, polygon )
394
395        assert allclose( res, [0,1,2,4,3] )
396        assert count == 3
397
398
399        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
400        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
401        res, count = separate_points_by_polygon( points, polygon )
402
403        assert allclose( res, [1,2,3,5,4,0] )       
404        assert count == 3
405       
406
407    def test_populate_polygon(self):
408
409        polygon = [[0,0], [1,0], [1,1], [0,1]]
410        points = populate_polygon(polygon, 5)
411
412        assert len(points) == 5
413        for point in points:
414            assert is_inside_polygon(point, polygon)
415
416
417        #Very convoluted polygon
418        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
419
420        points = populate_polygon(polygon, 5)
421
422        assert len(points) == 5
423        for point in points:
424            assert is_inside_polygon(point, polygon)
425
426
427    def test_populate_polygon_with_exclude(self):
428       
429
430        polygon = [[0,0], [1,0], [1,1], [0,1]]
431        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
432        points = populate_polygon(polygon, 5, exclude = [ex_poly])
433
434        assert len(points) == 5
435        for point in points:
436            assert is_inside_polygon(point, polygon)
437            assert not is_inside_polygon(point, ex_poly)           
438
439
440        #overlap
441        polygon = [[0,0], [1,0], [1,1], [0,1]]
442        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
443        points = populate_polygon(polygon, 5, exclude = [ex_poly])
444
445        assert len(points) == 5
446        for point in points:
447            assert is_inside_polygon(point, polygon)
448            assert not is_inside_polygon(point, ex_poly)                       
449       
450        #Multiple
451        polygon = [[0,0], [1,0], [1,1], [0,1]]
452        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
453        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter       
454       
455        points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])
456
457        assert len(points) == 20
458        for point in points:
459            assert is_inside_polygon(point, polygon)
460            assert not is_inside_polygon(point, ex_poly1)
461            assert not is_inside_polygon(point, ex_poly2)                               
462       
463
464        #Very convoluted polygon
465        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
466        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
467        points = populate_polygon(polygon, 20, exclude = [ex_poly])
468       
469        assert len(points) == 20
470        for point in points:
471            assert is_inside_polygon(point, polygon)
472            assert not is_inside_polygon(point, ex_poly), '%s' %str(point)                       
473
474
475    def test_populate_polygon_with_exclude2(self):
476       
477
478        min_outer = 0 
479        max_outer = 1000
480        polygon_outer = [[min_outer,min_outer],[max_outer,min_outer],
481                   [max_outer,max_outer],[min_outer,max_outer]]
482
483        delta = 10
484        min_inner1 = min_outer + delta
485        max_inner1 = max_outer - delta
486        inner1_polygon = [[min_inner1,min_inner1],[max_inner1,min_inner1],
487                   [max_inner1,max_inner1],[min_inner1,max_inner1]]
488     
489       
490        density_inner2 = 1000 
491        min_inner2 = min_outer +  2*delta
492        max_inner2 = max_outer -  2*delta
493        inner2_polygon = [[min_inner2,min_inner2],[max_inner2,min_inner2],
494                   [max_inner2,max_inner2],[min_inner2,max_inner2]]     
495       
496        points = populate_polygon(polygon_outer, 20, exclude = [inner1_polygon, inner2_polygon])
497
498        assert len(points) == 20
499        for point in points:
500            assert is_inside_polygon(point, polygon_outer)
501            assert not is_inside_polygon(point, inner1_polygon)
502            assert not is_inside_polygon(point, inner2_polygon)                               
503       
504
505        #Very convoluted polygon
506        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
507        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
508        points = populate_polygon(polygon, 20, exclude = [ex_poly])
509       
510        assert len(points) == 20
511        for point in points:
512            assert is_inside_polygon(point, polygon)
513            assert not is_inside_polygon(point, ex_poly), '%s' %str(point)                       
514
515    def test_point_in_polygon(self):
516        polygon = [[0,0], [1,0], [1,1], [0,1]]
517        point = point_in_polygon(polygon)
518        assert is_inside_polygon(point, polygon)
519
520        #this may get into a vicious loop
521        #polygon = [[1e32,1e54], [1,0], [1,1], [0,1]]
522       
523        polygon = [[1e15,1e7], [1,0], [1,1], [0,1]]
524        point = point_in_polygon(polygon)
525        assert is_inside_polygon(point, polygon)
526
527
528        polygon = [[0,0], [1,0], [1,1], [1e8,1e8]]
529        point = point_in_polygon(polygon)
530        assert is_inside_polygon(point, polygon)
531
532       
533        polygon = [[1e32,1e54], [-1e32,1e54], [1e32,-1e54]]
534        point = point_in_polygon(polygon)
535        assert is_inside_polygon(point, polygon)
536
537       
538        polygon = [[1e18,1e15], [1,0], [0,1]]
539        point = point_in_polygon(polygon)
540        assert is_inside_polygon(point, polygon)
541
542    def test_in_and_outside_polygon_main(self):
543
544
545        #Simplest case: Polygon is the unit square
546        polygon = [[0,0], [1,0], [1,1], [0,1]]
547
548        inside, outside =  in_and_outside_polygon( (0.5, 0.5), polygon )
549        assert inside[0] == 0
550        assert len(outside) == 0
551       
552        inside, outside =  in_and_outside_polygon(  (1., 0.5), polygon, closed=True)
553        assert inside[0] == 0
554        assert len(outside) == 0
555       
556        inside, outside =  in_and_outside_polygon(  (1., 0.5), polygon, closed=False)
557        assert len(inside) == 0
558        assert outside[0] == 0
559
560        points =  [(1., 0.25),(1., 0.75) ]
561        inside, outside =  in_and_outside_polygon( points, polygon, closed=True)
562        assert (inside, [0,1])
563        assert len(outside) == 0
564       
565        inside, outside =  in_and_outside_polygon( points, polygon, closed=False)
566        assert len(inside) == 0
567        assert (outside, [0,1])
568
569       
570        points =  [(100., 0.25),(0.5, 0.5) ] 
571        inside, outside =  in_and_outside_polygon( points, polygon)
572        assert (inside, [1])
573        assert outside[0] == 0
574       
575        points =  [(100., 0.25),(0.5, 0.5), (39,20), (0.6,0.7),(56,43),(67,90) ] 
576        inside, outside =  in_and_outside_polygon( points, polygon)
577        assert (inside, [1,3])
578        assert (outside, [0,2,4,5])
579
580
581    def test_intersection1(self):
582        line0 = [[-1,0], [1,0]]
583        line1 = [[0,-1], [0,1]]
584
585        p = intersection(line0, line1)
586        assert allclose(p, [0.0, 0.0])
587
588    def test_intersection2(self):
589        line0 = [[0,0], [24,12]]
590        line1 = [[0,12], [24,0]]
591
592        p = intersection(line0, line1)
593        assert allclose(p, [12.0, 6.0])
594
595        # Swap direction of one line
596        line1 = [[24,0], [0,12]]
597
598        p = intersection(line0, line1)
599        assert allclose(p, [12.0, 6.0])
600
601        # Swap order of lines
602        p = intersection(line1, line0)
603        assert allclose(p, [12.0, 6.0])       
604       
605    def test_intersection3(self):
606        line0 = [[0,0], [24,12]]
607        line1 = [[0,17], [24,0]]
608
609        p = intersection(line0, line1)
610        #print p
611        assert allclose(p, [14.068965517, 7.0344827586])
612
613        # Swap direction of one line
614        line1 = [[24,0], [0,17]]
615
616        p = intersection(line0, line1)
617        #print p
618        assert allclose(p, [14.068965517, 7.0344827586])       
619
620        # Swap order of lines
621        p = intersection(line1, line0)
622        assert allclose(p, [14.068965517, 7.0344827586])       
623
624
625    def test_intersection4(self):
626        line0 = [[0,0], [24,12]]
627        line1 = [[0,22], [21,0]]
628
629        p = intersection(line0, line1)
630        #print 'P',p
631
632       
633    def test_intersection_direction_invariance(self):
634        """This runs through a number of examples and checks that direction of lines don't matter.
635        """
636
637             
638        line0 = [[0,0], [100,100]]
639
640        common_end_point = [20, 150]
641       
642        for i in range(100):
643            x = 20 + i * 1.0/100
644
645            line1 = [[x,0], common_end_point]
646            p1 = intersection(line0, line1)
647
648            # Swap direction of line1
649            line1 = [common_end_point, [x,0]]           
650            p2 = intersection(line0, line1)
651
652            msg = 'Orientation of line shouldn not matter.\n'
653            msg += 'However, segment [%f,%f], [%f, %f]' %(x,
654                                                          0,
655                                                          common_end_point[0],
656                                                          common_end_point[1])
657            msg += ' gave %s, \nbut when reversed we got %s' %(p1, p2)
658            assert allclose(p1, p2), msg
659
660            # Swap order of lines
661            p3 = intersection(line1, line0)
662            msg = 'Order of lines gave different results'
663            assert allclose(p1, p3), msg
664
665    def test_no_intersection(self):
666        line0 = [[-1,1], [1,1]]
667        line1 = [[0,-1], [0,0]]
668
669        p = intersection(line0, line1)
670        assert p is None
671
672    def test_intersection_parallel(self):
673        line0 = [[-1,1], [1,1]]
674        line1 = [[-1,0], [5,0]]
675
676        p = intersection(line0, line1)
677        assert p is None               
678
679       
680       
681    def zzztest_inside_polygon_main(self):  \
682
683        #FIXME (Ole): Why is this disabled?
684        print "inside",inside
685        print "outside",outside
686       
687        assert not inside_polygon( (0.5, 1.5), polygon )
688        assert not inside_polygon( (0.5, -0.5), polygon )
689        assert not inside_polygon( (-0.5, 0.5), polygon )
690        assert not inside_polygon( (1.5, 0.5), polygon )
691
692        #Try point on borders
693        assert inside_polygon( (1., 0.5), polygon, closed=True)
694        assert inside_polygon( (0.5, 1), polygon, closed=True)
695        assert inside_polygon( (0., 0.5), polygon, closed=True)
696        assert inside_polygon( (0.5, 0.), polygon, closed=True)
697
698        assert not inside_polygon( (0.5, 1), polygon, closed=False)
699        assert not inside_polygon( (0., 0.5), polygon, closed=False)
700        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
701        assert not inside_polygon( (1., 0.5), polygon, closed=False)
702
703
704
705        #From real example (that failed)
706        polygon = [[20,20], [40,20], [40,40], [20,40]]
707        points = [ [40, 50] ]
708        res = inside_polygon(points, polygon)
709        assert len(res) == 0
710
711        polygon = [[20,20], [40,20], [40,40], [20,40]]
712        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
713        res = inside_polygon(points, polygon)
714        assert len(res) == 2
715        assert allclose(res, [0,1])
716
717    def test_polygon_area(self):
718
719        #Simplest case: Polygon is the unit square
720        polygon = [[0,0], [1,0], [1,1], [0,1]]
721        assert polygon_area(polygon) == 1
722
723        #Simple case: Polygon is a rectangle
724        polygon = [[0,0], [1,0], [1,4], [0,4]]
725        assert polygon_area(polygon) == 4
726
727        #Simple case: Polygon is a unit triangle
728        polygon = [[0,0], [1,0], [0,1]]
729        assert polygon_area(polygon) == 0.5
730
731        #Simple case: Polygon is a diamond
732        polygon = [[0,0], [1,1], [2,0], [1, -1]]
733        assert polygon_area(polygon) == 2.0
734
735    def test_poly_xy(self):
736 
737        #Simplest case: Polygon is the unit square
738        polygon = [[0,0], [1,0], [1,1], [0,1]]
739        x, y = poly_xy(polygon)
740        assert len(x) == len(polygon)+1
741        assert len(y) == len(polygon)+1
742        assert x[0] == 0
743        assert x[1] == 1
744        assert x[2] == 1
745        assert x[3] == 0
746        assert y[0] == 0
747        assert y[1] == 0
748        assert y[2] == 1
749        assert y[3] == 1
750
751        #Arbitrary polygon
752        polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]]
753        x, y = poly_xy(polygon)
754        assert len(x) == len(polygon)+1
755        assert len(y) == len(polygon)+1
756        assert x[0] == 1
757        assert x[1] == 1
758        assert x[2] == 100
759        assert x[3] == 1
760        assert x[4] == 3
761        assert y[0] == 5
762        assert y[1] == 1
763        assert y[2] == 10
764        assert y[3] == 10
765        assert y[4] == 6
766
767    # Disabled   
768    def xtest_plot_polygons(self):
769       
770        import os
771       
772        #Simplest case: Polygon is the unit square
773        polygon1 = [[0,0], [1,0], [1,1], [0,1]]
774        polygon2 = [[1,1], [2,1], [3,2], [2,2]]
775        v = plot_polygons([polygon1, polygon2],'test1')
776        assert len(v) == 4
777        assert v[0] == 0
778        assert v[1] == 3
779        assert v[2] == 0
780        assert v[3] == 2
781
782        #Another case
783        polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
784        v = plot_polygons([polygon2,polygon3],'test2')
785        assert len(v) == 4
786        assert v[0] == 1
787        assert v[1] == 100
788        assert v[2] == 1
789        assert v[3] == 10
790
791        os.remove('test1.png')
792        os.remove('test2.png')
793
794       
795    def test_inside_polygon_geospatial(self):
796
797
798        polygon_absolute = [[0,0], [1,0], [1,1], [0,1]]
799        poly_geo_ref = Geo_reference(57,100,100)
800       
801
802
803
804        #Simplest case: Polygon is the unit square
805        polygon_absolute = [[0,0], [1,0], [1,1], [0,1]]
806        poly_geo_ref = Geo_reference(57,100,100)
807        polygon = poly_geo_ref.change_points_geo_ref(polygon_absolute)
808        poly_spatial = Geospatial_data(polygon,
809                                       geo_reference=poly_geo_ref)
810       
811        points_absolute = (0.5, 0.5)
812        points_geo_ref = Geo_reference(57,78,-56)
813        points = points_geo_ref.change_points_geo_ref(points_absolute)
814        points_spatial = Geospatial_data(points,
815                                         geo_reference=points_geo_ref) 
816       
817        assert is_inside_polygon(points_absolute, polygon_absolute)
818        assert is_inside_polygon(ensure_numeric(points_absolute),
819                                 ensure_numeric(polygon_absolute))
820        assert is_inside_polygon(points_absolute, poly_spatial)
821        assert is_inside_polygon(points_spatial, poly_spatial)
822        assert is_inside_polygon(points_spatial, polygon_absolute)
823
824        assert is_inside_polygon(points_absolute, polygon_absolute)
825
826
827    def NOtest_decimate_polygon(self):
828
829        polygon = [[0,0], [10,10], [15,5], [20, 10],
830                   [25,0], [30,10], [40,-10], [35, -5]]
831
832        #plot_polygons([polygon], figname='test')
833       
834        dpoly = decimate_polygon(polygon, factor=2)
835
836        print dpoly
837       
838        assert len(dpoly)*2==len(polygon)
839
840        minx = maxx = polygon[0][0]
841        miny = maxy = polygon[0][1]
842        for point in polygon[1:]:
843            x, y = point
844           
845            if x < minx: minx = x
846            if x > maxx: maxx = x
847            if y < miny: miny = y
848            if y > maxy: maxy = y
849           
850
851        assert [minx, miny] in polygon
852        print minx, maxy
853        assert [minx, maxy] in polygon
854        assert [maxx, miny] in polygon
855        assert [maxx, maxy] in polygon               
856       
857
858       
859#-------------------------------------------------------------
860if __name__ == "__main__":
861    suite = unittest.makeSuite(Test_Polygon,'test')
862    #suite = unittest.makeSuite(Test_Polygon,'test_inside_polygon_geo_ref')
863    runner = unittest.TextTestRunner()
864    runner.run(suite)
865
866
867
868
Note: See TracBrowser for help on using the repository browser.