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

Last change on this file since 4424 was 4424, checked in by sexton, 17 years ago

incorporating new test for Polygon_function - used in ICs

File size: 23.1 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
8
9from polygon import *
10from anuga.coordinate_transforms.geo_reference import Geo_reference
11from anuga.geospatial_data.geospatial_data import Geospatial_data
12
13def test_function(x, y):
14    return x+y
15
16class Test_Polygon(unittest.TestCase):
17    def setUp(self):
18        pass
19
20    def tearDown(self):
21        pass
22
23
24    def test_that_C_extension_compiles(self):
25        FN = 'polygon_ext.c'
26        try:
27            import polygon_ext
28        except:
29            from compile import compile
30
31            try:
32                compile(FN)
33            except:
34                raise 'Could not compile %s' %FN
35            else:
36                import polygon_ext
37
38
39    #Polygon stuff
40    def test_polygon_function_constants(self):
41        p1 = [[0,0], [10,0], [10,10], [0,10]]
42        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
43
44        f = Polygon_function( [(p1, 1.0)] )
45        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
46        assert allclose(z, [1,1,0,0])
47
48
49        f = Polygon_function( [(p2, 2.0)] )
50        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
51        assert allclose(z, [2,0,0,2])
52
53
54        #Combined
55        f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
56        z = f([5, 5, 27, 35], [5, 9, 8, -5])
57        assert allclose(z, [2,1,0,2])
58
59    def test_polygon_function_csvfile(self):
60        from os import sep, getenv
61        #home = getenv('ANUGAHOME')
62        #p1 = read_polygon(home+sep+'anuga_core'+sep+'source'+sep+'anuga'+sep+'utilities'+sep+'mainland_only.csv')
63        home = getenv('INUNDATIONHOME')
64        p1 = read_polygon(home+sep+'data'+sep+'western_australia'+sep+'dampier_tsunami_scenario_2006'+sep+'anuga'+sep+'polygons'+sep+'2007polys'+sep+'mainland_only.csv')
65        #p1 = read_polygon('mainland_only.csv')
66       
67        f = Polygon_function( [(p1, 10.0)] )
68        z = f([430000,480000], [7720000, 7690000]) #first outside, second inside
69
70        assert allclose(z, [0,10])
71
72    def test_polygon_function_georef(self):
73        """Check that georeferencing works
74        """
75
76        from anuga.coordinate_transforms.geo_reference import Geo_reference
77
78        geo = Geo_reference(56, 200, 1000)
79
80        #Make points 'absolute'
81        p1 = [[200,1000], [210,1000], [210,1010], [200,1010]]
82        p2 = [[200,1000], [210,1010], [215,1005], [220, 1010], [225,1000],
83              [230,1010], [240,990]]
84
85        f = Polygon_function( [(p1, 1.0)], geo_reference = geo )
86        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
87
88        assert allclose(z, [1,1,0,0])
89
90
91        f = Polygon_function( [(p2, 2.0)], geo_reference = geo )
92        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
93        assert allclose(z, [2,0,0,2])
94
95
96        #Combined
97        f = Polygon_function( [(p1, 1.0), (p2, 2.0)], geo_reference = geo )
98        z = f([5, 5, 27, 35], [5, 9, 8, -5])
99        assert allclose(z, [2,1,0,2])
100
101
102        #Check that it would fail without geo
103        f = Polygon_function( [(p1, 1.0), (p2, 2.0)])
104        z = f([5, 5, 27, 35], [5, 9, 8, -5])
105        assert not allclose(z, [2,1,0,2])       
106
107
108
109    def test_polygon_function_callable(self):
110        """Check that values passed into Polygon_function can be callable
111        themselves.
112        """
113        p1 = [[0,0], [10,0], [10,10], [0,10]]
114        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
115
116        f = Polygon_function( [(p1, test_function)] )
117        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
118        assert allclose(z, [10,14,0,0])
119
120        #Combined
121        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
122        z = f([5, 5, 27, 35], [5, 9, 8, -5])
123        assert allclose(z, [2,14,0,2])
124
125
126        #Combined w default
127        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
128        z = f([5, 5, 27, 35], [5, 9, 8, -5])
129        assert allclose(z, [2,14,3.14,2])
130
131
132        #Combined w default func
133        f = Polygon_function( [(p1, test_function), (p2, 2.0)],
134                              default = test_function)
135        z = f([5, 5, 27, 35], [5, 9, 8, -5])
136        assert allclose(z, [2,14,35,2])
137
138
139
140    def test_point_on_line(self):
141
142        #Endpoints first
143        assert point_on_line( 0, 0, 0,0, 1,0 )
144        assert point_on_line( 1, 0, 0,0, 1,0 )
145
146        #Then points on line
147        assert point_on_line( 0.5, 0, 0,0, 1,0 )
148        assert point_on_line( 0, 0.5, 0,1, 0,0 )
149        assert point_on_line( 1, 0.5, 1,1, 1,0 )
150        assert point_on_line( 0.5, 0.5, 0,0, 1,1 )
151
152        #Then points not on line
153        assert not point_on_line( 0.5, 0, 0,1, 1,1 )
154        assert not point_on_line( 0, 0.5, 0,0, 1,1 )
155
156        #From real example that failed
157        assert not point_on_line( 40,50, 40,20, 40,40 )
158
159
160        #From real example that failed
161        assert not point_on_line( 40,19, 40,20, 40,40 )
162
163
164
165
166    def test_is_inside_polygon_main(self):
167
168
169        #Simplest case: Polygon is the unit square
170        polygon = [[0,0], [1,0], [1,1], [0,1]]
171
172        assert is_inside_polygon( (0.5, 0.5), polygon )
173        assert not is_inside_polygon( (0.5, 1.5), polygon )
174        assert not is_inside_polygon( (0.5, -0.5), polygon )
175        assert not is_inside_polygon( (-0.5, 0.5), polygon )
176        assert not is_inside_polygon( (1.5, 0.5), polygon )
177
178        #Try point on borders
179        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
180        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
181        assert is_inside_polygon( (0., 0.5), polygon, closed=True)
182        assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
183
184        assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
185        assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
186        assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
187        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
188
189
190    def test_inside_polygon_main(self):
191
192        #Simplest case: Polygon is the unit square
193        polygon = [[0,0], [1,0], [1,1], [0,1]]       
194
195        #From real example (that failed)
196        polygon = [[20,20], [40,20], [40,40], [20,40]]
197        points = [ [40, 50] ]
198        res = inside_polygon(points, polygon)
199        assert len(res) == 0
200
201        polygon = [[20,20], [40,20], [40,40], [20,40]]
202        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
203        res = inside_polygon(points, polygon)
204        assert len(res) == 2
205        assert allclose(res, [0,1])
206
207
208
209        #More convoluted and non convex polygon
210        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
211        assert is_inside_polygon( (0.5, 0.5), polygon )
212        assert is_inside_polygon( (1, -0.5), polygon )
213        assert is_inside_polygon( (1.5, 0), polygon )
214
215        assert not is_inside_polygon( (0.5, 1.5), polygon )
216        assert not is_inside_polygon( (0.5, -0.5), polygon )
217
218
219        #Very convoluted polygon
220        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
221        assert is_inside_polygon( (5, 5), polygon )
222        assert is_inside_polygon( (17, 7), polygon )
223        assert is_inside_polygon( (27, 2), polygon )
224        assert is_inside_polygon( (35, -5), polygon )
225        assert not is_inside_polygon( (15, 7), polygon )
226        assert not is_inside_polygon( (24, 3), polygon )
227        assert not is_inside_polygon( (25, -10), polygon )
228
229
230
231        #Another combination (that failed)
232        polygon = [[0,0], [10,0], [10,10], [0,10]]
233        assert is_inside_polygon( (5, 5), polygon )
234        assert is_inside_polygon( (7, 7), polygon )
235        assert not is_inside_polygon( (-17, 7), polygon )
236        assert not is_inside_polygon( (7, 17), polygon )
237        assert not is_inside_polygon( (17, 7), polygon )
238        assert not is_inside_polygon( (27, 8), polygon )
239        assert not is_inside_polygon( (35, -5), polygon )
240
241
242
243
244        #Multiple polygons
245
246        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
247                   [10,10], [11,10], [11,11], [10,11], [10,10]]
248        assert is_inside_polygon( (0.5, 0.5), polygon )
249        assert is_inside_polygon( (10.5, 10.5), polygon )
250
251        #FIXME: Fails if point is 5.5, 5.5
252        assert not is_inside_polygon( (0, 5.5), polygon )
253
254        #Polygon with a hole
255        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
256                   [0,0], [1,0], [1,1], [0,1], [0,0]]
257
258        assert is_inside_polygon( (0, -0.5), polygon )
259        assert not is_inside_polygon( (0.5, 0.5), polygon )
260
261
262
263    def test_duplicate_points_being_ok(self):
264
265
266        #Simplest case: Polygon is the unit square
267        polygon = [[0,0], [1,0], [1,0], [1,0], [1,1], [0,1], [0,0]]
268
269        assert is_inside_polygon( (0.5, 0.5), polygon )
270        assert not is_inside_polygon( (0.5, 1.5), polygon )
271        assert not is_inside_polygon( (0.5, -0.5), polygon )
272        assert not is_inside_polygon( (-0.5, 0.5), polygon )
273        assert not is_inside_polygon( (1.5, 0.5), polygon )
274
275        #Try point on borders
276        assert is_inside_polygon( (1., 0.5), polygon, closed=True)
277        assert is_inside_polygon( (0.5, 1), polygon, closed=True)
278        assert is_inside_polygon( (0., 0.5), polygon, closed=True)
279        assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
280
281        assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
282        assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
283        assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
284        assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
285
286        #From real example
287        polygon = [[20,20], [40,20], [40,40], [20,40]]
288        points = [ [40, 50] ]
289        res = inside_polygon(points, polygon)
290        assert len(res) == 0
291
292       
293
294    def test_inside_polygon_vector_version(self):
295        #Now try the vector formulation returning indices
296        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
297        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
298        res = inside_polygon( points, polygon, verbose=False )
299
300        assert allclose( res, [0,1,2] )
301
302    def test_outside_polygon(self):
303        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
304
305        assert not is_outside_polygon( [0.5, 0.5], U )
306        #evaluate to False as the point 0.5, 0.5 is inside the unit square
307       
308        assert is_outside_polygon( [1.5, 0.5], U )
309        #evaluate to True as the point 1.5, 0.5 is outside the unit square
310       
311        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
312        assert allclose( indices, [1] )
313       
314        #One more test of vector formulation returning indices
315        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
316        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
317        res = outside_polygon( points, polygon )
318
319        assert allclose( res, [3, 4] )
320
321
322
323        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
324        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
325        res = outside_polygon( points, polygon )
326
327        assert allclose( res, [0, 4, 5] )       
328     
329    def test_outside_polygon2(self):
330        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
331   
332        assert not outside_polygon( [0.5, 1.0], U, closed = True )
333        #evaluate to False as the point 0.5, 1.0 is inside the unit square
334       
335        assert is_outside_polygon( [0.5, 1.0], U, closed = False )
336        #evaluate to True as the point 0.5, 1.0 is outside the unit square
337
338    def test_all_outside_polygon(self):
339        """Test case where all points are outside poly
340        """
341       
342        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
343
344        points = [[2,2], [1,3], [-1,1], [0,2]] #All outside
345
346
347        indices, count = separate_points_by_polygon(points, U)
348        #print indices, count
349        assert count == 0 #None inside
350        assert allclose(indices, [3,2,1,0])
351
352        indices = outside_polygon(points, U, closed = True)
353        assert allclose(indices, [0,1,2,3])
354
355        indices = inside_polygon(points, U, closed = True)
356        assert allclose(indices, [])               
357
358
359    def test_all_inside_polygon(self):
360        """Test case where all points are inside poly
361        """
362       
363        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
364
365        points = [[0.5,0.5], [0.2,0.3], [0,0.5]] #All inside (or on edge)
366
367
368        indices, count = separate_points_by_polygon(points, U)
369        assert count == 3 #All inside
370        assert allclose(indices, [0,1,2])
371
372        indices = outside_polygon(points, U, closed = True)
373        assert allclose(indices, [])
374
375        indices = inside_polygon(points, U, closed = True)
376        assert allclose(indices, [0,1,2])
377       
378
379    def test_separate_points_by_polygon(self):
380        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
381
382        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
383        assert allclose( indices, [0,2,1] )
384        assert count == 2
385       
386        #One more test of vector formulation returning indices
387        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
388        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
389        res, count = separate_points_by_polygon( points, polygon )
390
391        assert allclose( res, [0,1,2,4,3] )
392        assert count == 3
393
394
395        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
396        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
397        res, count = separate_points_by_polygon( points, polygon )
398
399        assert allclose( res, [1,2,3,5,4,0] )       
400        assert count == 3
401       
402
403    def test_populate_polygon(self):
404
405        polygon = [[0,0], [1,0], [1,1], [0,1]]
406        points = populate_polygon(polygon, 5)
407
408        assert len(points) == 5
409        for point in points:
410            assert is_inside_polygon(point, polygon)
411
412
413        #Very convoluted polygon
414        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
415
416        points = populate_polygon(polygon, 5)
417
418        assert len(points) == 5
419        for point in points:
420            assert is_inside_polygon(point, polygon)
421
422
423    def test_populate_polygon_with_exclude(self):
424       
425
426        polygon = [[0,0], [1,0], [1,1], [0,1]]
427        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
428        points = populate_polygon(polygon, 5, exclude = [ex_poly])
429
430        assert len(points) == 5
431        for point in points:
432            assert is_inside_polygon(point, polygon)
433            assert not is_inside_polygon(point, ex_poly)           
434
435
436        #overlap
437        polygon = [[0,0], [1,0], [1,1], [0,1]]
438        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
439        points = populate_polygon(polygon, 5, exclude = [ex_poly])
440
441        assert len(points) == 5
442        for point in points:
443            assert is_inside_polygon(point, polygon)
444            assert not is_inside_polygon(point, ex_poly)                       
445       
446        #Multiple
447        polygon = [[0,0], [1,0], [1,1], [0,1]]
448        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
449        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter       
450       
451        points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])
452
453        assert len(points) == 20
454        for point in points:
455            assert is_inside_polygon(point, polygon)
456            assert not is_inside_polygon(point, ex_poly1)
457            assert not is_inside_polygon(point, ex_poly2)                               
458       
459
460        #Very convoluted polygon
461        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
462        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
463        points = populate_polygon(polygon, 20, exclude = [ex_poly])
464       
465        assert len(points) == 20
466        for point in points:
467            assert is_inside_polygon(point, polygon)
468            assert not is_inside_polygon(point, ex_poly), '%s' %str(point)                       
469
470
471    def test_populate_polygon_with_exclude2(self):
472       
473
474        min_outer = 0 
475        max_outer = 1000
476        polygon_outer = [[min_outer,min_outer],[max_outer,min_outer],
477                   [max_outer,max_outer],[min_outer,max_outer]]
478
479        delta = 10
480        min_inner1 = min_outer + delta
481        max_inner1 = max_outer - delta
482        inner1_polygon = [[min_inner1,min_inner1],[max_inner1,min_inner1],
483                   [max_inner1,max_inner1],[min_inner1,max_inner1]]
484     
485       
486        density_inner2 = 1000 
487        min_inner2 = min_outer +  2*delta
488        max_inner2 = max_outer -  2*delta
489        inner2_polygon = [[min_inner2,min_inner2],[max_inner2,min_inner2],
490                   [max_inner2,max_inner2],[min_inner2,max_inner2]]     
491       
492        points = populate_polygon(polygon_outer, 20, exclude = [inner1_polygon, inner2_polygon])
493
494        assert len(points) == 20
495        for point in points:
496            assert is_inside_polygon(point, polygon_outer)
497            assert not is_inside_polygon(point, inner1_polygon)
498            assert not is_inside_polygon(point, inner2_polygon)                               
499       
500
501        #Very convoluted polygon
502        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
503        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
504        points = populate_polygon(polygon, 20, exclude = [ex_poly])
505       
506        assert len(points) == 20
507        for point in points:
508            assert is_inside_polygon(point, polygon)
509            assert not is_inside_polygon(point, ex_poly), '%s' %str(point)                       
510
511    def test_point_in_polygon(self):
512        polygon = [[0,0], [1,0], [1,1], [0,1]]
513        point = point_in_polygon(polygon)
514        assert is_inside_polygon(point, polygon)
515
516        #this may get into a vicious loop
517        #polygon = [[1e32,1e54], [1,0], [1,1], [0,1]]
518       
519        polygon = [[1e15,1e7], [1,0], [1,1], [0,1]]
520        point = point_in_polygon(polygon)
521        assert is_inside_polygon(point, polygon)
522
523
524        polygon = [[0,0], [1,0], [1,1], [1e8,1e8]]
525        point = point_in_polygon(polygon)
526        assert is_inside_polygon(point, polygon)
527
528       
529        polygon = [[1e32,1e54], [-1e32,1e54], [1e32,-1e54]]
530        point = point_in_polygon(polygon)
531        assert is_inside_polygon(point, polygon)
532
533       
534        polygon = [[1e18,1e15], [1,0], [0,1]]
535        point = point_in_polygon(polygon)
536        assert is_inside_polygon(point, polygon)
537
538    def test_in_and_outside_polygon_main(self):
539
540
541        #Simplest case: Polygon is the unit square
542        polygon = [[0,0], [1,0], [1,1], [0,1]]
543
544        inside, outside =  in_and_outside_polygon( (0.5, 0.5), polygon )
545        assert inside[0] == 0
546        assert len(outside) == 0
547       
548        inside, outside =  in_and_outside_polygon(  (1., 0.5), polygon, closed=True)
549        assert inside[0] == 0
550        assert len(outside) == 0
551       
552        inside, outside =  in_and_outside_polygon(  (1., 0.5), polygon, closed=False)
553        assert len(inside) == 0
554        assert outside[0] == 0
555
556        points =  [(1., 0.25),(1., 0.75) ]
557        inside, outside =  in_and_outside_polygon( points, polygon, closed=True)
558        assert (inside, [0,1])
559        assert len(outside) == 0
560       
561        inside, outside =  in_and_outside_polygon( points, polygon, closed=False)
562        assert len(inside) == 0
563        assert (outside, [0,1])
564
565       
566        points =  [(100., 0.25),(0.5, 0.5) ] 
567        inside, outside =  in_and_outside_polygon( points, polygon)
568        assert (inside, [1])
569        assert outside[0] == 0
570       
571        points =  [(100., 0.25),(0.5, 0.5), (39,20), (0.6,0.7),(56,43),(67,90) ] 
572        inside, outside =  in_and_outside_polygon( points, polygon)
573        assert (inside, [1,3])
574        assert (outside, [0,2,4,5])
575       
576    def zzztest_inside_polygon_main(self): 
577        print "inside",inside
578        print "outside",outside
579       
580        assert not inside_polygon( (0.5, 1.5), polygon )
581        assert not inside_polygon( (0.5, -0.5), polygon )
582        assert not inside_polygon( (-0.5, 0.5), polygon )
583        assert not inside_polygon( (1.5, 0.5), polygon )
584
585        #Try point on borders
586        assert inside_polygon( (1., 0.5), polygon, closed=True)
587        assert inside_polygon( (0.5, 1), polygon, closed=True)
588        assert inside_polygon( (0., 0.5), polygon, closed=True)
589        assert inside_polygon( (0.5, 0.), polygon, closed=True)
590
591        assert not inside_polygon( (0.5, 1), polygon, closed=False)
592        assert not inside_polygon( (0., 0.5), polygon, closed=False)
593        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
594        assert not inside_polygon( (1., 0.5), polygon, closed=False)
595
596
597
598        #From real example (that failed)
599        polygon = [[20,20], [40,20], [40,40], [20,40]]
600        points = [ [40, 50] ]
601        res = inside_polygon(points, polygon)
602        assert len(res) == 0
603
604        polygon = [[20,20], [40,20], [40,40], [20,40]]
605        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
606        res = inside_polygon(points, polygon)
607        assert len(res) == 2
608        assert allclose(res, [0,1])
609
610    def test_polygon_area(self):
611
612        #Simplest case: Polygon is the unit square
613        polygon = [[0,0], [1,0], [1,1], [0,1]]
614        assert polygon_area(polygon) == 1
615
616        #Simple case: Polygon is a rectangle
617        polygon = [[0,0], [1,0], [1,4], [0,4]]
618        assert polygon_area(polygon) == 4
619
620        #Simple case: Polygon is a unit triangle
621        polygon = [[0,0], [1,0], [0,1]]
622        assert polygon_area(polygon) == 0.5
623
624        #Simple case: Polygon is a diamond
625        polygon = [[0,0], [1,1], [2,0], [1, -1]]
626        assert polygon_area(polygon) == 2.0
627
628    def test_poly_xy(self):
629 
630        #Simplest case: Polygon is the unit square
631        polygon = [[0,0], [1,0], [1,1], [0,1]]
632        x, y = poly_xy(polygon)
633        assert len(x) == len(polygon)+1
634        assert len(y) == len(polygon)+1
635        assert x[0] == 0
636        assert x[1] == 1
637        assert x[2] == 1
638        assert x[3] == 0
639        assert y[0] == 0
640        assert y[1] == 0
641        assert y[2] == 1
642        assert y[3] == 1
643
644        #Arbitrary polygon
645        polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]]
646        x, y = poly_xy(polygon)
647        assert len(x) == len(polygon)+1
648        assert len(y) == len(polygon)+1
649        assert x[0] == 1
650        assert x[1] == 1
651        assert x[2] == 100
652        assert x[3] == 1
653        assert x[4] == 3
654        assert y[0] == 5
655        assert y[1] == 1
656        assert y[2] == 10
657        assert y[3] == 10
658        assert y[4] == 6
659
660    # Disabled   
661    def xtest_plot_polygons(self):
662       
663        import os
664       
665        #Simplest case: Polygon is the unit square
666        polygon1 = [[0,0], [1,0], [1,1], [0,1]]
667        polygon2 = [[1,1], [2,1], [3,2], [2,2]]
668        v = plot_polygons([polygon1, polygon2],'test1')
669        assert len(v) == 4
670        assert v[0] == 0
671        assert v[1] == 3
672        assert v[2] == 0
673        assert v[3] == 2
674
675        #Another case
676        polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
677        v = plot_polygons([polygon2,polygon3],'test2')
678        assert len(v) == 4
679        assert v[0] == 1
680        assert v[1] == 100
681        assert v[2] == 1
682        assert v[3] == 10
683
684        os.remove('test1.png')
685        os.remove('test2.png')
686
687       
688    def test_inside_polygon_geospatial(self):
689
690
691        polygon_absolute = [[0,0], [1,0], [1,1], [0,1]]
692        poly_geo_ref = Geo_reference(57,100,100)
693       
694
695
696
697        #Simplest case: Polygon is the unit square
698        polygon_absolute = [[0,0], [1,0], [1,1], [0,1]]
699        poly_geo_ref = Geo_reference(57,100,100)
700        polygon = poly_geo_ref.change_points_geo_ref(polygon_absolute)
701        poly_spatial = Geospatial_data(polygon,
702                                       geo_reference=poly_geo_ref)
703       
704        points_absolute = (0.5, 0.5)
705        points_geo_ref = Geo_reference(57,78,-56)
706        points = points_geo_ref.change_points_geo_ref(points_absolute)
707        points_spatial = Geospatial_data(points,
708                                         geo_reference=points_geo_ref) 
709       
710        assert is_inside_polygon(points_absolute, polygon_absolute)
711        assert is_inside_polygon(ensure_numeric(points_absolute),
712                                 ensure_numeric(polygon_absolute))
713        assert is_inside_polygon(points_absolute, poly_spatial)
714        assert is_inside_polygon(points_spatial, poly_spatial)
715        assert is_inside_polygon(points_spatial, polygon_absolute)
716
717        assert is_inside_polygon(points_absolute, polygon_absolute)
718       
719       
720#-------------------------------------------------------------
721if __name__ == "__main__":
722    suite = unittest.makeSuite(Test_Polygon,'test')
723    #suite = unittest.makeSuite(Test_Polygon,'test_inside_polygon_geo_ref')
724    runner = unittest.TextTestRunner()
725    runner.run(suite)
726
727
728
729
Note: See TracBrowser for help on using the repository browser.