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

Last change on this file since 4589 was 4528, checked in by ole, 18 years ago

Fixed polygon test so that it can run both from within its own directory and from test_all.py one level above it.

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