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

Last change on this file since 5079 was 5072, checked in by ole, 17 years ago

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