source: inundation/utilities/test_polygon.py @ 2655

Last change on this file since 2655 was 2655, checked in by duncan, 18 years ago

Adding ability to remove points outside of the mesh when interpolating

File size: 17.2 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
6from math import sqrt, pi
7
8from polygon import *
9
10
11def test_function(x, y):
12    return x+y
13
14class Test_Polygon(unittest.TestCase):
15    def setUp(self):
16        pass
17
18    def tearDown(self):
19        pass
20
21
22    def test_that_C_extension_compiles(self):
23        FN = 'polygon_ext.c'
24        try:
25            import polygon_ext
26        except:
27            from compile import compile
28
29            try:
30                compile(FN)
31            except:
32                raise 'Could not compile %s' %FN
33            else:
34                import polygon_ext
35
36
37    #Polygon stuff
38    def test_polygon_function_constants(self):
39        p1 = [[0,0], [10,0], [10,10], [0,10]]
40        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
41
42        f = Polygon_function( [(p1, 1.0)] )
43        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
44        assert allclose(z, [1,1,0,0])
45
46
47        f = Polygon_function( [(p2, 2.0)] )
48        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
49        assert allclose(z, [2,0,0,2])
50
51
52        #Combined
53        f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
54        z = f([5, 5, 27, 35], [5, 9, 8, -5])
55        assert allclose(z, [2,1,0,2])
56
57
58    def test_polygon_function_georef(self):
59        """Check that georeferencing works
60        """
61
62        from coordinate_transforms.geo_reference import Geo_reference
63
64        geo = Geo_reference(56, 200, 1000)
65
66        #Make points 'absolute'
67        p1 = [[200,1000], [210,1000], [210,1010], [200,1010]]
68        p2 = [[200,1000], [210,1010], [215,1005], [220, 1010], [225,1000],
69              [230,1010], [240,990]]
70
71        f = Polygon_function( [(p1, 1.0)], geo_reference = geo )
72        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
73
74        assert allclose(z, [1,1,0,0])
75
76
77        f = Polygon_function( [(p2, 2.0)], geo_reference = geo )
78        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
79        assert allclose(z, [2,0,0,2])
80
81
82        #Combined
83        f = Polygon_function( [(p1, 1.0), (p2, 2.0)], geo_reference = geo )
84        z = f([5, 5, 27, 35], [5, 9, 8, -5])
85        assert allclose(z, [2,1,0,2])
86
87
88        #Check that it would fail without geo
89        f = Polygon_function( [(p1, 1.0), (p2, 2.0)])
90        z = f([5, 5, 27, 35], [5, 9, 8, -5])
91        assert not allclose(z, [2,1,0,2])       
92
93
94
95    def test_polygon_function_callable(self):
96        """Check that values passed into Polygon_function can be callable
97        themselves.
98        """
99        p1 = [[0,0], [10,0], [10,10], [0,10]]
100        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
101
102        f = Polygon_function( [(p1, test_function)] )
103        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
104        assert allclose(z, [10,14,0,0])
105
106        #Combined
107        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
108        z = f([5, 5, 27, 35], [5, 9, 8, -5])
109        assert allclose(z, [2,14,0,2])
110
111
112        #Combined w default
113        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
114        z = f([5, 5, 27, 35], [5, 9, 8, -5])
115        assert allclose(z, [2,14,3.14,2])
116
117
118        #Combined w default func
119        f = Polygon_function( [(p1, test_function), (p2, 2.0)],
120                              default = test_function)
121        z = f([5, 5, 27, 35], [5, 9, 8, -5])
122        assert allclose(z, [2,14,35,2])
123
124
125
126    def test_point_on_line(self):
127
128        #Endpoints first
129        assert point_on_line( 0, 0, 0,0, 1,0 )
130        assert point_on_line( 1, 0, 0,0, 1,0 )
131
132        #Then points on line
133        assert point_on_line( 0.5, 0, 0,0, 1,0 )
134        assert point_on_line( 0, 0.5, 0,1, 0,0 )
135        assert point_on_line( 1, 0.5, 1,1, 1,0 )
136        assert point_on_line( 0.5, 0.5, 0,0, 1,1 )
137
138        #Then points not on line
139        assert not point_on_line( 0.5, 0, 0,1, 1,1 )
140        assert not point_on_line( 0, 0.5, 0,0, 1,1 )
141
142        #From real example that failed
143        assert not point_on_line( 40,50, 40,20, 40,40 )
144
145
146        #From real example that failed
147        assert not point_on_line( 40,19, 40,20, 40,40 )
148
149
150
151
152    def test_inside_polygon_main(self):
153
154
155        #Simplest case: Polygon is the unit square
156        polygon = [[0,0], [1,0], [1,1], [0,1]]
157
158        assert inside_polygon( (0.5, 0.5), polygon )
159        assert not inside_polygon( (0.5, 1.5), polygon )
160        assert not inside_polygon( (0.5, -0.5), polygon )
161        assert not inside_polygon( (-0.5, 0.5), polygon )
162        assert not inside_polygon( (1.5, 0.5), polygon )
163
164        #Try point on borders
165        assert inside_polygon( (1., 0.5), polygon, closed=True)
166        assert inside_polygon( (0.5, 1), polygon, closed=True)
167        assert inside_polygon( (0., 0.5), polygon, closed=True)
168        assert inside_polygon( (0.5, 0.), polygon, closed=True)
169
170        assert not inside_polygon( (0.5, 1), polygon, closed=False)
171        assert not inside_polygon( (0., 0.5), polygon, closed=False)
172        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
173        assert not inside_polygon( (1., 0.5), polygon, closed=False)
174
175
176
177        #From real example (that failed)
178        polygon = [[20,20], [40,20], [40,40], [20,40]]
179        points = [ [40, 50] ]
180        res = inside_polygon(points, polygon)
181        assert len(res) == 0
182
183        polygon = [[20,20], [40,20], [40,40], [20,40]]
184        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
185        res = inside_polygon(points, polygon)
186        assert len(res) == 2
187        assert allclose(res, [0,1])
188
189
190
191        #More convoluted and non convex polygon
192        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
193        assert inside_polygon( (0.5, 0.5), polygon )
194        assert inside_polygon( (1, -0.5), polygon )
195        assert inside_polygon( (1.5, 0), polygon )
196
197        assert not inside_polygon( (0.5, 1.5), polygon )
198        assert not inside_polygon( (0.5, -0.5), polygon )
199
200
201        #Very convoluted polygon
202        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
203        assert inside_polygon( (5, 5), polygon )
204        assert inside_polygon( (17, 7), polygon )
205        assert inside_polygon( (27, 2), polygon )
206        assert inside_polygon( (35, -5), polygon )
207        assert not inside_polygon( (15, 7), polygon )
208        assert not inside_polygon( (24, 3), polygon )
209        assert not inside_polygon( (25, -10), polygon )
210
211
212
213        #Another combination (that failed)
214        polygon = [[0,0], [10,0], [10,10], [0,10]]
215        assert inside_polygon( (5, 5), polygon )
216        assert inside_polygon( (7, 7), polygon )
217        assert not inside_polygon( (-17, 7), polygon )
218        assert not inside_polygon( (7, 17), polygon )
219        assert not inside_polygon( (17, 7), polygon )
220        assert not inside_polygon( (27, 8), polygon )
221        assert not inside_polygon( (35, -5), polygon )
222
223
224
225
226        #Multiple polygons
227
228        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
229                   [10,10], [11,10], [11,11], [10,11], [10,10]]
230        assert inside_polygon( (0.5, 0.5), polygon )
231        assert inside_polygon( (10.5, 10.5), polygon )
232
233        #FIXME: Fails if point is 5.5, 5.5
234        assert not inside_polygon( (0, 5.5), polygon )
235
236        #Polygon with a hole
237        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
238                   [0,0], [1,0], [1,1], [0,1], [0,0]]
239
240        assert inside_polygon( (0, -0.5), polygon )
241        assert not inside_polygon( (0.5, 0.5), polygon )
242
243
244
245    def test_duplicate_points_being_ok(self):
246
247
248        #Simplest case: Polygon is the unit square
249        polygon = [[0,0], [1,0], [1,0], [1,0], [1,1], [0,1], [0,0]]
250
251        assert inside_polygon( (0.5, 0.5), polygon )
252        assert not inside_polygon( (0.5, 1.5), polygon )
253        assert not inside_polygon( (0.5, -0.5), polygon )
254        assert not inside_polygon( (-0.5, 0.5), polygon )
255        assert not inside_polygon( (1.5, 0.5), polygon )
256
257        #Try point on borders
258        assert inside_polygon( (1., 0.5), polygon, closed=True)
259        assert inside_polygon( (0.5, 1), polygon, closed=True)
260        assert inside_polygon( (0., 0.5), polygon, closed=True)
261        assert inside_polygon( (0.5, 0.), polygon, closed=True)
262
263        assert not inside_polygon( (0.5, 1), polygon, closed=False)
264        assert not inside_polygon( (0., 0.5), polygon, closed=False)
265        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
266        assert not inside_polygon( (1., 0.5), polygon, closed=False)
267
268        #From real example
269        polygon = [[20,20], [40,20], [40,40], [20,40]]
270        points = [ [40, 50] ]
271        res = inside_polygon(points, polygon)
272        assert len(res) == 0
273
274       
275
276    def test_inside_polygon_vector_version(self):
277        #Now try the vector formulation returning indices
278        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
279        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
280        res = inside_polygon( points, polygon, verbose=False )
281
282        assert allclose( res, [0,1,2] )
283
284    def test_outside_polygon(self):
285        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
286
287        assert not outside_polygon( [0.5, 0.5], U )
288        #evaluate to False as the point 0.5, 0.5 is inside the unit square
289       
290        assert outside_polygon( [1.5, 0.5], U )
291        #evaluate to True as the point 1.5, 0.5 is outside the unit square
292       
293        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
294        assert allclose( indices, [1] )
295       
296        #One more test of vector formulation returning indices
297        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
298        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
299        res = outside_polygon( points, polygon )
300
301        assert allclose( res, [3, 4] )
302
303
304
305        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
306        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
307        res = outside_polygon( points, polygon )
308
309        assert allclose( res, [0, 4, 5] )       
310     
311    def test_outside_polygon2(self):
312        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
313   
314        assert not outside_polygon( [0.5, 1.0], U, closed = True )
315        #evaluate to False as the point 0.5, 1.0 is inside the unit square
316       
317        assert outside_polygon( [0.5, 1.0], U, closed = False )
318        #evaluate to True as the point 0.5, 1.0 is outside the unit square
319
320    def test_all_outside_polygon(self):
321        """Test case where all points are outside poly
322        """
323       
324        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
325
326        points = [[2,2], [1,3], [-1,1], [0,2]] #All outside
327
328
329        indices, count = separate_points_by_polygon(points, U)
330        #print indices, count
331        assert count == 0 #None inside
332        assert allclose(indices, [3,2,1,0])
333
334        indices = outside_polygon(points, U, closed = True)
335        assert allclose(indices, [0,1,2,3])
336
337        indices = inside_polygon(points, U, closed = True)
338        assert allclose(indices, [])               
339
340
341    def test_all_inside_polygon(self):
342        """Test case where all points are inside poly
343        """
344       
345        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
346
347        points = [[0.5,0.5], [0.2,0.3], [0,0.5]] #All inside (or on edge)
348
349
350        indices, count = separate_points_by_polygon(points, U)
351        assert count == 3 #All inside
352        assert allclose(indices, [0,1,2])
353
354        indices = outside_polygon(points, U, closed = True)
355        assert allclose(indices, [])
356
357        indices = inside_polygon(points, U, closed = True)
358        assert allclose(indices, [0,1,2])
359       
360
361    def test_separate_points_by_polygon(self):
362        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
363
364        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
365        assert allclose( indices, [0,2,1] )
366        assert count == 2
367       
368        #One more test of vector formulation returning indices
369        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
370        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
371        res, count = separate_points_by_polygon( points, polygon )
372
373        assert allclose( res, [0,1,2,4,3] )
374        assert count == 3
375
376
377        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
378        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
379        res, count = separate_points_by_polygon( points, polygon )
380
381        assert allclose( res, [1,2,3,5,4,0] )       
382        assert count == 3
383       
384
385    def test_populate_polygon(self):
386
387        polygon = [[0,0], [1,0], [1,1], [0,1]]
388        points = populate_polygon(polygon, 5)
389
390        assert len(points) == 5
391        for point in points:
392            assert inside_polygon(point, polygon)
393
394
395        #Very convoluted polygon
396        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
397
398        points = populate_polygon(polygon, 5)
399
400        assert len(points) == 5
401        for point in points:
402            assert inside_polygon(point, polygon)
403
404
405    def test_populate_polygon_with_exclude(self):
406       
407
408        polygon = [[0,0], [1,0], [1,1], [0,1]]
409        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
410        points = populate_polygon(polygon, 5, exclude = [ex_poly])
411
412        assert len(points) == 5
413        for point in points:
414            assert inside_polygon(point, polygon)
415            assert not inside_polygon(point, ex_poly)           
416
417
418        #overlap
419        polygon = [[0,0], [1,0], [1,1], [0,1]]
420        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
421        points = populate_polygon(polygon, 5, exclude = [ex_poly])
422
423        assert len(points) == 5
424        for point in points:
425            assert inside_polygon(point, polygon)
426            assert not inside_polygon(point, ex_poly)                       
427       
428        #Multiple
429        polygon = [[0,0], [1,0], [1,1], [0,1]]
430        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
431        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter       
432       
433        points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])
434
435        assert len(points) == 20
436        for point in points:
437            assert inside_polygon(point, polygon)
438            assert not inside_polygon(point, ex_poly1)
439            assert not inside_polygon(point, ex_poly2)                               
440       
441
442        #Very convoluted polygon
443        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
444        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
445        points = populate_polygon(polygon, 20, exclude = [ex_poly])
446       
447        assert len(points) == 20
448        for point in points:
449            assert inside_polygon(point, polygon)
450            assert not inside_polygon(point, ex_poly), '%s' %str(point)                       
451
452    def test_point_in_polygon(self):
453        polygon = [[0,0], [1,0], [1,1], [0,1]]
454        point = point_in_polygon(polygon)
455        assert inside_polygon(point, polygon)
456
457        #this may get into a vicious loop
458        #polygon = [[1e32,1e54], [1,0], [1,1], [0,1]]
459       
460        polygon = [[1e15,1e7], [1,0], [1,1], [0,1]]
461        point = point_in_polygon(polygon)
462        assert inside_polygon(point, polygon)
463
464
465        polygon = [[0,0], [1,0], [1,1], [1e8,1e8]]
466        point = point_in_polygon(polygon)
467        assert inside_polygon(point, polygon)
468
469       
470        polygon = [[1e32,1e54], [-1e32,1e54], [1e32,-1e54]]
471        point = point_in_polygon(polygon)
472        assert inside_polygon(point, polygon)
473
474       
475        polygon = [[1e18,1e15], [1,0], [0,1]]
476        point = point_in_polygon(polygon)
477        assert inside_polygon(point, polygon)
478
479    def test_in_and_outside_polygon_main(self):
480
481
482        #Simplest case: Polygon is the unit square
483        polygon = [[0,0], [1,0], [1,1], [0,1]]
484
485        inside, outside =  in_and_outside_polygon( (0.5, 0.5), polygon )
486        assert inside[0] == 0
487        assert len(outside) == 0
488       
489        inside, outside =  in_and_outside_polygon(  (1., 0.5), polygon, closed=True)
490        assert inside[0] == 0
491        assert len(outside) == 0
492       
493        inside, outside =  in_and_outside_polygon(  (1., 0.5), polygon, closed=False)
494        assert len(inside) == 0
495        assert outside[0] == 0
496
497        points =  [(1., 0.25),(1., 0.75) ]
498        inside, outside =  in_and_outside_polygon( points, polygon, closed=True)
499        assert (inside, [0,1])
500        assert len(outside) == 0
501       
502        inside, outside =  in_and_outside_polygon( points, polygon, closed=False)
503        assert len(inside) == 0
504        assert (outside, [0,1])
505
506       
507        points =  [(100., 0.25),(0.5, 0.5) ] 
508        inside, outside =  in_and_outside_polygon( points, polygon)
509        assert (inside, [1])
510        assert outside[0] == 0
511       
512        points =  [(100., 0.25),(0.5, 0.5), (39,20), (0.6,0.7),(56,43),(67,90) ] 
513        inside, outside =  in_and_outside_polygon( points, polygon)
514        assert (inside, [1,3])
515        assert (outside, [0,2,4,5])
516       
517    def zzztest_inside_polygon_main(self): 
518        print "inside",inside
519        print "outside",outside
520       
521        assert not inside_polygon( (0.5, 1.5), polygon )
522        assert not inside_polygon( (0.5, -0.5), polygon )
523        assert not inside_polygon( (-0.5, 0.5), polygon )
524        assert not inside_polygon( (1.5, 0.5), polygon )
525
526        #Try point on borders
527        assert inside_polygon( (1., 0.5), polygon, closed=True)
528        assert inside_polygon( (0.5, 1), polygon, closed=True)
529        assert inside_polygon( (0., 0.5), polygon, closed=True)
530        assert inside_polygon( (0.5, 0.), polygon, closed=True)
531
532        assert not inside_polygon( (0.5, 1), polygon, closed=False)
533        assert not inside_polygon( (0., 0.5), polygon, closed=False)
534        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
535        assert not inside_polygon( (1., 0.5), polygon, closed=False)
536
537
538
539        #From real example (that failed)
540        polygon = [[20,20], [40,20], [40,40], [20,40]]
541        points = [ [40, 50] ]
542        res = inside_polygon(points, polygon)
543        assert len(res) == 0
544
545        polygon = [[20,20], [40,20], [40,40], [20,40]]
546        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
547        res = inside_polygon(points, polygon)
548        assert len(res) == 2
549        assert allclose(res, [0,1])
550       
551#-------------------------------------------------------------
552if __name__ == "__main__":
553    suite = unittest.makeSuite(Test_Polygon,'test')
554    runner = unittest.TextTestRunner()
555    runner.run(suite)
556
557
558
559
Note: See TracBrowser for help on using the repository browser.