source: inundation/utilities/test_polygon.py @ 2378

Last change on this file since 2378 was 2375, checked in by ole, 19 years ago

Added georeferencing to polygon_function and added test

File size: 13.7 KB
RevLine 
[1910]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
[2375]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
[1910]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
[2375]125
[1910]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    def test_inside_polygon_vector_version(self):
244        #Now try the vector formulation returning indices
245        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
246        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
247        res = inside_polygon( points, polygon, verbose=False )
248
249        assert allclose( res, [0,1,2] )
250
251    def test_outside_polygon(self):
252        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
253
254        assert not outside_polygon( [0.5, 0.5], U )
255        #evaluate to False as the point 0.5, 0.5 is inside the unit square
256       
257        assert outside_polygon( [1.5, 0.5], U )
258        #evaluate to True as the point 1.5, 0.5 is outside the unit square
259       
260        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
261        assert allclose( indices, [1] )
262       
263        #One more test of vector formulation returning indices
264        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
265        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
266        res = outside_polygon( points, polygon )
267
268        assert allclose( res, [3, 4] )
269
270
271
272        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
273        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
274        res = outside_polygon( points, polygon )
275
276        assert allclose( res, [0, 4, 5] )       
277     
278    def test_outside_polygon2(self):
279        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
280   
281        assert not outside_polygon( [0.5, 1.0], U, closed = True )
282        #evaluate to False as the point 0.5, 1.0 is inside the unit square
283       
284        assert outside_polygon( [0.5, 1.0], U, closed = False )
285        #evaluate to True as the point 0.5, 1.0 is outside the unit square
286
[2062]287    def test_all_outside_polygon(self):
288        """Test case where all points are outside poly
289        """
290       
291        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
292
293        points = [[2,2], [1,3], [-1,1], [0,2]] #All outside
294
295
296        indices, count = separate_points_by_polygon(points, U)
297        #print indices, count
298        assert count == 0 #None inside
299        assert allclose(indices, [3,2,1,0])
300
301        indices = outside_polygon(points, U, closed = True)
302        assert allclose(indices, [0,1,2,3])
303
304        indices = inside_polygon(points, U, closed = True)
305        assert allclose(indices, [])               
306
307
308    def test_all_inside_polygon(self):
309        """Test case where all points are inside poly
310        """
311       
312        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
313
314        points = [[0.5,0.5], [0.2,0.3], [0,0.5]] #All inside (or on edge)
315
316
317        indices, count = separate_points_by_polygon(points, U)
318        assert count == 3 #All inside
319        assert allclose(indices, [0,1,2])
320
321        indices = outside_polygon(points, U, closed = True)
322        assert allclose(indices, [])
323
324        indices = inside_polygon(points, U, closed = True)
325        assert allclose(indices, [0,1,2])
326       
327
[1910]328    def test_separate_points_by_polygon(self):
329        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
330
331        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
332        assert allclose( indices, [0,2,1] )
333        assert count == 2
334       
335        #One more test of vector formulation returning indices
336        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
337        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
338        res, count = separate_points_by_polygon( points, polygon )
339
340        assert allclose( res, [0,1,2,4,3] )
341        assert count == 3
342
343
344        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
345        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
346        res, count = separate_points_by_polygon( points, polygon )
347
348        assert allclose( res, [1,2,3,5,4,0] )       
349        assert count == 3
350       
351
352    def test_populate_polygon(self):
353
354        polygon = [[0,0], [1,0], [1,1], [0,1]]
355        points = populate_polygon(polygon, 5)
356
357        assert len(points) == 5
358        for point in points:
359            assert inside_polygon(point, polygon)
360
361
362        #Very convoluted polygon
363        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
364
365        points = populate_polygon(polygon, 5)
366
367        assert len(points) == 5
368        for point in points:
369            assert inside_polygon(point, polygon)
370
371
372    def test_populate_polygon_with_exclude(self):
373       
374
375        polygon = [[0,0], [1,0], [1,1], [0,1]]
376        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
377        points = populate_polygon(polygon, 5, exclude = [ex_poly])
378
379        assert len(points) == 5
380        for point in points:
381            assert inside_polygon(point, polygon)
382            assert not inside_polygon(point, ex_poly)           
383
384
385        #overlap
386        polygon = [[0,0], [1,0], [1,1], [0,1]]
387        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
388        points = populate_polygon(polygon, 5, exclude = [ex_poly])
389
390        assert len(points) == 5
391        for point in points:
392            assert inside_polygon(point, polygon)
393            assert not inside_polygon(point, ex_poly)                       
394       
395        #Multiple
396        polygon = [[0,0], [1,0], [1,1], [0,1]]
397        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
398        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter       
399       
400        points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])
401
402        assert len(points) == 20
403        for point in points:
404            assert inside_polygon(point, polygon)
405            assert not inside_polygon(point, ex_poly1)
406            assert not inside_polygon(point, ex_poly2)                               
407       
408
409        #Very convoluted polygon
410        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
411        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
412        points = populate_polygon(polygon, 20, exclude = [ex_poly])
413       
414        assert len(points) == 20
415        for point in points:
416            assert inside_polygon(point, polygon)
417            assert not inside_polygon(point, ex_poly), '%s' %str(point)                       
418
[2200]419    def test_point_in_polygon(self):
420        polygon = [[0,0], [1,0], [1,1], [0,1]]
421        point = point_in_polygon(polygon)
422        assert inside_polygon(point, polygon)
[1910]423
[2310]424        #this may get into a vicious loop
425        #polygon = [[1e32,1e54], [1,0], [1,1], [0,1]]
426       
[2311]427        polygon = [[1e15,1e7], [1,0], [1,1], [0,1]]
[2200]428        point = point_in_polygon(polygon)
429        assert inside_polygon(point, polygon)
[1910]430
[2200]431
432        polygon = [[0,0], [1,0], [1,1], [1e8,1e8]]
433        point = point_in_polygon(polygon)
434        assert inside_polygon(point, polygon)
[2311]435
[2200]436       
[2310]437        polygon = [[1e32,1e54], [-1e32,1e54], [1e32,-1e54]]
[2200]438        point = point_in_polygon(polygon)
439        assert inside_polygon(point, polygon)
[2311]440
[2200]441       
[2311]442        polygon = [[1e18,1e15], [1,0], [0,1]]
[2200]443        point = point_in_polygon(polygon)
444        assert inside_polygon(point, polygon)
445
[1910]446#-------------------------------------------------------------
447if __name__ == "__main__":
448    suite = unittest.makeSuite(Test_Polygon,'test')
449    runner = unittest.TextTestRunner()
450    runner.run(suite)
451
452
453
454
Note: See TracBrowser for help on using the repository browser.