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

Last change on this file since 5940 was 5940, checked in by rwilson, 15 years ago

Modify the Rajaraman test case.

File size: 43.5 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        f = Polygon_function( [(p1, 10.0)] )
72        z = f([430000,480000], [490000,7720000]) # first outside, second inside
73       
74        assert allclose(z, [0,10])
75
76    def test_polygon_function_georef(self):
77        """Check that georeferencing works
78        """
79
80        from anuga.coordinate_transforms.geo_reference import Geo_reference
81
82        geo = Geo_reference(56, 200, 1000)
83
84        # Make points 'absolute'
85        p1 = [[200,1000], [210,1000], [210,1010], [200,1010]]
86        p2 = [[200,1000], [210,1010], [215,1005], [220, 1010], [225,1000],
87              [230,1010], [240,990]]
88
89        f = Polygon_function( [(p1, 1.0)], geo_reference=geo)
90        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
91
92        assert allclose(z, [1,1,0,0])
93
94
95        f = Polygon_function( [(p2, 2.0)], geo_reference=geo)
96        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
97        assert allclose(z, [2,0,0,2])
98
99
100        # Combined
101        f = Polygon_function( [(p1, 1.0), (p2, 2.0)], geo_reference=geo)
102        z = f([5, 5, 27, 35], [5, 9, 8, -5])
103        assert allclose(z, [2,1,0,2])
104
105
106        # Check that it would fail without geo
107        f = Polygon_function( [(p1, 1.0), (p2, 2.0)])
108        z = f([5, 5, 27, 35], [5, 9, 8, -5])
109        assert not allclose(z, [2,1,0,2])       
110
111
112
113    def test_polygon_function_callable(self):
114        """Check that values passed into Polygon_function can be callable
115        themselves.
116        """
117        p1 = [[0,0], [10,0], [10,10], [0,10]]
118        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
119
120        f = Polygon_function( [(p1, test_function)] )
121        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
122        assert allclose(z, [10,14,0,0])
123
124        # Combined
125        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
126        z = f([5, 5, 27, 35], [5, 9, 8, -5])
127        assert allclose(z, [2,14,0,2])
128
129
130        # Combined w default
131        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
132        z = f([5, 5, 27, 35], [5, 9, 8, -5])
133        assert allclose(z, [2,14,3.14,2])
134
135
136        # Combined w default func
137        f = Polygon_function( [(p1, test_function), (p2, 2.0)],
138                              default = test_function)
139        z = f([5, 5, 27, 35], [5, 9, 8, -5])
140        assert allclose(z, [2,14,35,2])
141
142
143
144    def test_point_on_line(self):
145
146        # Endpoints first
147        assert point_on_line( [0, 0], [[0,0], [1,0]] )
148        assert point_on_line( [1, 0], [[0,0], [1,0]] )
149
150        # Then points on line
151        assert point_on_line( [0.5, 0], [[0,0], [1,0]] )
152        assert point_on_line( [0, 0.5], [[0,1], [0,0]] )
153        assert point_on_line( [1, 0.5], [[1,1], [1,0]] )
154        assert point_on_line( [0.5, 0.5], [[0,0], [1,1]] )
155
156        # Then points not on line
157        assert not point_on_line( [0.5, 0], [[0,1], [1,1]] )
158        assert not point_on_line( [0, 0.5], [[0,0], [1,1]] )
159
160        # From real example that failed
161        assert not point_on_line( [40,50], [[40,20], [40,40]] )
162
163
164        # From real example that failed
165        assert not point_on_line( [40,19], [[40,20], [40,40]] )
166
167        # Degenerate line
168        assert point_on_line( [40,19], [[40,19], [40,19]] )
169        assert not point_on_line( [40,19], [[40,40], [40,40]] )       
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
584    def test_intersection1(self):
585        line0 = [[-1,0], [1,0]]
586        line1 = [[0,-1], [0,1]]
587
588        status, value = intersection(line0, line1)
589        assert status == 1
590        assert allclose(value, [0.0, 0.0])
591
592    def test_intersection2(self):
593        line0 = [[0,0], [24,12]]
594        line1 = [[0,12], [24,0]]
595
596        status, value = intersection(line0, line1)
597        assert status == 1
598        assert allclose(value, [12.0, 6.0])
599
600        # Swap direction of one line
601        line1 = [[24,0], [0,12]]
602
603        status, value = intersection(line0, line1)
604        assert status == 1
605        assert allclose(value, [12.0, 6.0])
606
607        # Swap order of lines
608        status, value = intersection(line1, line0)
609        assert status == 1
610        assert allclose(value, [12.0, 6.0])       
611       
612    def test_intersection3(self):
613        line0 = [[0,0], [24,12]]
614        line1 = [[0,17], [24,0]]
615
616        status, value = intersection(line0, line1)
617        assert status == 1
618        assert allclose(value, [14.068965517, 7.0344827586])
619
620        # Swap direction of one line
621        line1 = [[24,0], [0,17]]
622
623        status, value = intersection(line0, line1)
624        assert status == 1
625        assert allclose(value, [14.068965517, 7.0344827586])       
626
627        # Swap order of lines
628        status, value = intersection(line1, line0)
629        assert status == 1       
630        assert allclose(value, [14.068965517, 7.0344827586])       
631
632
633    def test_intersection_endpoints(self):
634        """test_intersection_endpoints(self):
635
636        Test that coinciding endpoints are picked up
637        """
638        line0 = [[0,0], [1,1]]
639        line1 = [[1,1], [2,1]]
640
641        status, value = intersection(line0, line1)
642        assert status == 1
643        assert allclose(value, [1.0, 1.0])
644
645
646        line0 = [[1,1], [2,0]]
647        line1 = [[1,1], [2,1]]
648
649        status, value = intersection(line0, line1)
650        assert status == 1
651        assert allclose(value, [1.0, 1.0])       
652       
653
654    # This function is a helper function for the test_intersection_bug_20081110_?()
655    # set of tests.
656    # This function should never be run directly by the unittest code.
657    def helper_test_parallel_intersection_code(self, P1, P2, P3, P4):
658        # lines in same direction, no overlap
659        # 0:         ---->----
660        # 1:                     --------->-----------
661        line0 = [P1,P2]
662        line1 = [P3,P4]
663        status, value = intersection(line0, line1)
664        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
665                               (str(status), str(value)))
666        self.failUnless(value is None, 'Expected value of None, got %s' %
667                                       str(value))
668
669        # lines in same direction, no overlap
670        # 0:         ----<----
671        # 1:                     ---------<-----------
672        line0 = [P2,P1]
673        line1 = [P4,P3]
674        status, value = intersection(line0, line1)
675        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
676                               (str(status), str(value)))
677        self.failUnless(value is None, 'Expected value of None, got %s' %
678                                       str(value))
679
680        # lines in opposite direction, no overlap
681        # 0:         ----<----
682        # 1:                     --------->-----------
683        line0 = [P2,P1]
684        line1 = [P3,P4]
685        status, value = intersection(line0, line1)
686        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
687                               (str(status), str(value)))
688        self.failUnless(value is None, 'Expected value of None, got %s' %
689                                       str(value))
690
691        # lines in opposite direction, no overlap
692        # 0:         ---->----
693        # 1:                     ---------<-----------
694        line0 = [P1,P2]
695        line1 = [P4,P3]
696        status, value = intersection(line0, line1)
697        self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
698                               (str(status), str(value)))
699        self.failUnless(value is None, 'Expected value of None, got %s' %
700                                       str(value))
701
702        # line0 fully within line1, same direction
703        # 0:         ---->----
704        # 1:    --------->-----------
705        # value should be line0:
706        #            ---->----
707        line0 = [P2,P3]
708        line1 = [P1,P4]
709        status, value = intersection(line0, line1)
710        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
711                               (str(status), str(value)))
712        self.failUnless(allclose(value, line0))
713
714        # line0 fully within line1, same direction
715        # 0:         ----<----
716        # 1:    ---------<-----------
717        # value should be line0:
718        #            ----<----
719        line0 = [P3,P2]
720        line1 = [P4,P1]
721        status, value = intersection(line0, line1)
722        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
723                               (str(status), str(value)))
724        self.failUnless(allclose(value, line0))
725
726        # line0 fully within line1, opposite direction
727        # 0:         ----<----
728        # 1:    --------->-----------
729        # value should be line0:
730        #            ----<----
731        line0 = [P3,P2]
732        line1 = [P1,P4]
733        status, value = intersection(line0, line1)
734        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
735                               (str(status), str(value)))
736        self.failUnless(allclose(value, line0))
737
738        # line0 fully within line1, opposite direction
739        # 0:         ---->----
740        # 1:    ---------<-----------
741        # value should be line0:
742        #            ---->----
743        line0 = [P2,P3]
744        line1 = [P4,P1]
745        status, value = intersection(line0, line1)
746        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
747                               (str(status), str(value)))
748        self.failUnless(allclose(value, line0))
749
750        # line1 fully within line0, same direction
751        # 0:    --------->-----------
752        # 1:         ---->----
753        # value should be line1:
754        #            ---->----
755        line0 = [P1,P4]
756        line1 = [P2,P3]
757        status, value = intersection(line0, line1)
758        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
759                               (str(status), str(value)))
760        self.failUnless(allclose(value, line1))
761
762        # line1 fully within line0, same direction
763        # 0:    ---------<-----------
764        # 1:         ----<----
765        # value should be line1:
766        #            ----<----
767        line0 = [P4,P1]
768        line1 = [P3,P2]
769        status, value = intersection(line0, line1)
770        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
771                               (str(status), str(value)))
772        self.failUnless(allclose(value, line1))
773
774        # line1 fully within line0, opposite direction
775        # 0:    ---------<-----------
776        # 1:         ---->----
777        # value should be line1:
778        #            ---->----
779        line0 = [P4,P1]
780        line1 = [P2,P3]
781        status, value = intersection(line0, line1)
782        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
783                               (str(status), str(value)))
784        self.failUnless(allclose(value, line1))
785
786        # line1 fully within line0, opposite direction
787        # 0:    --------->-----------
788        # 1:         ----<----
789        # value should be line1:
790        #            ----<----
791        line0 = [P1,P4]
792        line1 = [P3,P2]
793        status, value = intersection(line0, line1)
794        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
795                               (str(status), str(value)))
796        self.failUnless(allclose(value, line1))
797
798        # line in same direction, partial overlap
799        # 0:    ----->-----
800        # 1:       ------->--------
801        # value should be segment line1_start->line0_end:
802        #          --->----
803        line0 = [P1,P3]
804        line1 = [P2,P4]
805        status, value = intersection(line0, line1)
806        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
807                               (str(status), str(value)))
808        self.failUnless(allclose(value, [line1[0],line0[1]]))
809
810        # line in same direction, partial overlap
811        # 0:    -----<-----
812        # 1:       -------<--------
813        # value should be segment line0_start->line1_end:
814        #          ---<----
815        line0 = [P3,P1]
816        line1 = [P4,P2]
817        status, value = intersection(line0, line1)
818        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
819                               (str(status), str(value)))
820        self.failUnless(allclose(value, [line0[0],line1[1]]))
821
822        # line in opposite direction, partial overlap
823        # 0:    -----<-----
824        # 1:       ------->--------
825        # value should be segment line0_start->line1_start:
826        #          ---<----
827        line0 = [P3,P1]
828        line1 = [P2,P4]
829        status, value = intersection(line0, line1)
830        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
831                               (str(status), str(value)))
832        self.failUnless(allclose(value, [line0[0],line1[0]]))
833
834        # line in opposite direction, partial overlap
835        # 0:    ----->-----
836        # 1:       -------<--------
837        # value should be segment line1_end->line0_end:
838        #          --->----
839        line0 = [P1,P3]
840        line1 = [P4,P2]
841        status, value = intersection(line0, line1)
842        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
843                               (str(status), str(value)))
844        self.failUnless(allclose(value, [line1[1],line0[1]]))
845
846        # line in same direction, partial overlap
847        # 0:       ------>------
848        # 1:    ------>------
849        # value should be segment line0_start->line1_end:
850        #          --->----
851        line0 = [P2,P4]
852        line1 = [P1,P3]
853        status, value = intersection(line0, line1)
854        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
855                               (str(status), str(value)))
856        self.failUnless(allclose(value, [line0[0],line1[1]]))
857
858        # line in same direction, partial overlap
859        # 0:       ------<------
860        # 1:    ------<------
861        # value should be segment line1_start->line0_end:
862        #          ----<-----
863        line0 = [P4,P2]
864        line1 = [P3,P1]
865        status, value = intersection(line0, line1)
866        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
867                               (str(status), str(value)))
868        self.failUnless(allclose(value, [line1[0],line0[1]]))
869
870        # line in opposite direction, partial overlap
871        # 0:       ------<------
872        # 1:    ----->------
873        # value should be segment line1_end->line0_end:
874        #          --->----
875        line0 = [P4,P2]
876        line1 = [P1,P3]
877        status, value = intersection(line0, line1)
878        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
879                               (str(status), str(value)))
880        self.failUnless(allclose(value, [line1[1],line0[1]]))
881
882        # line in opposite direction, partial overlap
883        # 0:       ------>------
884        # 1:    -----<------
885        # value should be segment line0_start->line1_start:
886        #          ---<----
887        line0 = [P2,P4]
888        line1 = [P3,P1]
889        status, value = intersection(line0, line1)
890        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
891                               (str(status), str(value)))
892        self.failUnless(allclose(value, [line0[0],line1[0]]))
893
894       
895    def test_intersection_bug_20081110_1(self):
896        """test_intersection_bug_20081110(self)
897
898        Run through all possible cases of parallel lines.
899        """
900
901        # define 4 collinear points
902        #    P1---P2---P3---P4
903        P1 = [1.0, 0.0]
904        P2 = [2.0, 0.0]
905        P3 = [3.0, 0.0]
906        P4 = [4.0, 0.0]
907       
908        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)       
909
910
911    def test_intersection_bug_20081110_2(self):
912        """test_intersection_bug_20081110(self)
913
914        Run through all possible cases of parallel lines.
915        """
916
917        # define 4 *almost* collinear points
918        #    P1---P2---P3---P4
919        P1 = [1.0, 1.0e-9]
920        P2 = [2.0, 0.0]
921        P3 = [3.0, 0.0]
922        P4 = [4.0, 0.0]
923
924        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)       
925
926
927    def test_intersection_bug_20081110_3(self):
928        """test_intersection_bug_20081110(self)
929
930        Run through all possible cases of parallel lines.
931        """
932
933        # define 4 *almost* collinear points
934        #    P1---P2---P3---P4
935        P1 = [1.0, 0.0]
936        P2 = [2.0, 1.0e-9]
937        P3 = [3.0, 0.0]
938        P4 = [4.0, 0.0]
939
940        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)       
941
942
943    def test_intersection_bug_20081110_4(self):
944        """test_intersection_bug_20081110(self)
945
946        Run through all possible cases of parallel lines.
947        """
948
949        # define 4 *almost* collinear points
950        #    P1---P2---P3---P4
951        P1 = [1.0, 0.0]
952        P2 = [2.0, 0.0]
953        P3 = [3.0, 1.0e-9]
954        P4 = [4.0, 0.0]
955
956        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)       
957
958
959    def test_intersection_bug_20081110_5(self):
960        """test_intersection_bug_20081110(self)
961
962        Run through all possible cases of parallel lines.
963        """
964
965        # define 4 *almost* collinear points
966        #    P1---P2---P3---P4
967        P1 = [1.0, 0.0]
968        P2 = [2.0, 0.0]
969        P3 = [3.0, 0.0]
970        P4 = [4.0, 1.0e-9]
971
972        self.helper_test_parallel_intersection_code(P1, P2, P3, P4)       
973
974
975    def test_intersection_bug_20081111(self):
976        """test_intersection_bug_20081111(self)
977
978        Case from Rajaraman that is failing still.
979        """
980
981        # define 4  collinear points
982        #
983        #   Y
984        #   ^
985        #   |         x (2,2)
986        #   |        /
987        #   |       /
988        #   |      /
989        #   |     /
990        #   |    x (1,1)
991        #   |   /
992        #   |  /
993        #   | /
994        #   |/
995        #   x------------>X
996
997        line0 = [[2.0, 2.0],[0.0,0.0]]
998        line1 = [[1.0, 1.0],[0.0,0.0]]
999
1000        status, value = intersection(line0, line1)
1001        print 'status=%s, value=%s' % (str(status), str(value))
1002       
1003        self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1004                               (str(status), str(value)))
1005        self.failUnless(allclose(value, [line0[0],line1[1]]))
1006
1007
1008    def test_intersection_direction_invariance(self):
1009        """This runs through a number of examples and checks that direction of lines don't matter.
1010        """
1011             
1012        line0 = [[0,0], [100,100]]
1013
1014        common_end_point = [20, 150]
1015       
1016        for i in range(100):
1017            x = 20 + i * 1.0/100
1018
1019            line1 = [[x,0], common_end_point]
1020            status, p1 = intersection(line0, line1)
1021            assert status == 1
1022
1023
1024            # Swap direction of line1
1025            line1 = [common_end_point, [x,0]]           
1026            status, p2 = intersection(line0, line1)
1027            assert status == 1           
1028
1029            msg = 'Orientation of line shouldn not matter.\n'
1030            msg += 'However, segment [%f,%f], [%f, %f]' %(x,
1031                                                          0,
1032                                                          common_end_point[0],
1033                                                          common_end_point[1])
1034            msg += ' gave %s, \nbut when reversed we got %s' %(p1, p2)
1035            assert allclose(p1, p2), msg
1036
1037            # Swap order of lines
1038            status, p3 = intersection(line1, line0)
1039            assert status == 1                       
1040            msg = 'Order of lines gave different results'
1041            assert allclose(p1, p3), msg
1042           
1043
1044    def test_no_intersection(self):
1045        line0 = [[-1,1], [1,1]]
1046        line1 = [[0,-1], [0,0]]
1047
1048        status, value = intersection(line0, line1)
1049        assert status == 0
1050        assert value is None
1051       
1052
1053    def test_intersection_parallel(self):
1054        line0 = [[-1,1], [1,1]]
1055        line1 = [[-1,0], [5,0]]
1056
1057        status, value = intersection(line0, line1)
1058        assert status == 4       
1059        assert value is None
1060
1061
1062        line0 = [[0,0], [10,100]]
1063        line1 = [[-10,5], [0,105]]
1064
1065        status, value = intersection(line0, line1)
1066        assert status == 4               
1067        assert value is None       
1068
1069
1070    def test_intersection_coincide(self):
1071        """def test_intersection_coincide(self):
1072        Test what happens whe two lines partly coincide
1073        """
1074
1075        # Overlap 1
1076        line0 = [[0,0], [5,0]]
1077        line1 = [[-3,0], [3,0]]
1078
1079        status, value = intersection(line0, line1)
1080        assert status == 2
1081        assert allclose(value, [[0,0], [3,0]])
1082
1083        # Overlap 2
1084        line0 = [[-10,0], [5,0]]
1085        line1 = [[-3,0], [10,0]]
1086
1087        status, value = intersection(line0, line1)
1088        assert status == 2
1089        assert allclose(value, [[-3, 0], [5,0]])       
1090
1091        # Inclusion 1
1092        line0 = [[0,0], [5,0]]
1093        line1 = [[2,0], [3,0]]
1094
1095        status, value = intersection(line0, line1)
1096        assert status == 2       
1097        assert allclose(value, line1)
1098
1099        # Inclusion 2
1100        line0 = [[1,0], [5,0]]
1101        line1 = [[-10,0], [15,0]]
1102
1103        status, value = intersection(line0, line1)
1104        assert status == 2       
1105        assert allclose(value, line0)                                       
1106
1107
1108        # Exclusion (no intersection)
1109        line0 = [[-10,0], [1,0]]
1110        line1 = [[3,0], [15,0]]
1111
1112        status, value = intersection(line0, line1)
1113        assert status == 3       
1114        assert value is None
1115       
1116
1117        # Try examples with some slope (y=2*x+5)
1118
1119        # Overlap
1120        line0 = [[0,5], [7,19]]
1121        line1 = [[1,7], [10,25]]
1122        status, value = intersection(line0, line1)
1123        assert status == 2               
1124        assert allclose(value, [[1, 7], [7, 19]])
1125
1126        status, value = intersection(line1, line0)
1127        assert status == 2
1128        assert allclose(value, [[1, 7], [7, 19]])
1129
1130        # Swap direction
1131        line0 = [[7,19], [0,5]]
1132        line1 = [[1,7], [10,25]]
1133        status, value = intersection(line0, line1)
1134        assert status == 2
1135        assert allclose(value, [[7, 19], [1, 7]])
1136
1137        line0 = [[0,5], [7,19]]
1138        line1 = [[10,25], [1,7]]
1139        status, value = intersection(line0, line1)
1140        assert status == 2
1141        assert allclose(value, [[1, 7], [7, 19]])       
1142       
1143
1144        # Inclusion
1145        line0 = [[1,7], [7,19]]
1146        line1 = [[0,5], [10,25]]
1147        status, value = intersection(line0, line1)
1148        assert status == 2                       
1149        assert allclose(value, [[1,7], [7, 19]])               
1150
1151        line0 = [[0,5], [10,25]]
1152        line1 = [[1,7], [7,19]]
1153        status, value = intersection(line0, line1)
1154        assert status == 2                       
1155        assert allclose(value, [[1,7], [7, 19]])
1156
1157
1158        line0 = [[0,5], [10,25]]
1159        line1 = [[7,19], [1,7]]
1160        status, value = intersection(line0, line1)
1161        assert status == 2                       
1162        assert allclose(value, [[7, 19], [1, 7]])                       
1163       
1164       
1165    def zzztest_inside_polygon_main(self):  \
1166
1167        #FIXME (Ole): Why is this disabled?
1168        print "inside",inside
1169        print "outside",outside
1170       
1171        assert not inside_polygon( (0.5, 1.5), polygon )
1172        assert not inside_polygon( (0.5, -0.5), polygon )
1173        assert not inside_polygon( (-0.5, 0.5), polygon )
1174        assert not inside_polygon( (1.5, 0.5), polygon )
1175
1176        #Try point on borders
1177        assert inside_polygon( (1., 0.5), polygon, closed=True)
1178        assert inside_polygon( (0.5, 1), polygon, closed=True)
1179        assert inside_polygon( (0., 0.5), polygon, closed=True)
1180        assert inside_polygon( (0.5, 0.), polygon, closed=True)
1181
1182        assert not inside_polygon( (0.5, 1), polygon, closed=False)
1183        assert not inside_polygon( (0., 0.5), polygon, closed=False)
1184        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
1185        assert not inside_polygon( (1., 0.5), polygon, closed=False)
1186
1187
1188
1189        #From real example (that failed)
1190        polygon = [[20,20], [40,20], [40,40], [20,40]]
1191        points = [ [40, 50] ]
1192        res = inside_polygon(points, polygon)
1193        assert len(res) == 0
1194
1195        polygon = [[20,20], [40,20], [40,40], [20,40]]
1196        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
1197        res = inside_polygon(points, polygon)
1198        assert len(res) == 2
1199        assert allclose(res, [0,1])
1200
1201    def test_polygon_area(self):
1202
1203        #Simplest case: Polygon is the unit square
1204        polygon = [[0,0], [1,0], [1,1], [0,1]]
1205        assert polygon_area(polygon) == 1
1206
1207        #Simple case: Polygon is a rectangle
1208        polygon = [[0,0], [1,0], [1,4], [0,4]]
1209        assert polygon_area(polygon) == 4
1210
1211        #Simple case: Polygon is a unit triangle
1212        polygon = [[0,0], [1,0], [0,1]]
1213        assert polygon_area(polygon) == 0.5
1214
1215        #Simple case: Polygon is a diamond
1216        polygon = [[0,0], [1,1], [2,0], [1, -1]]
1217        assert polygon_area(polygon) == 2.0
1218
1219    def test_poly_xy(self):
1220 
1221        #Simplest case: Polygon is the unit square
1222        polygon = [[0,0], [1,0], [1,1], [0,1]]
1223        x, y = poly_xy(polygon)
1224        assert len(x) == len(polygon)+1
1225        assert len(y) == len(polygon)+1
1226        assert x[0] == 0
1227        assert x[1] == 1
1228        assert x[2] == 1
1229        assert x[3] == 0
1230        assert y[0] == 0
1231        assert y[1] == 0
1232        assert y[2] == 1
1233        assert y[3] == 1
1234
1235        #Arbitrary polygon
1236        polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]]
1237        x, y = poly_xy(polygon)
1238        assert len(x) == len(polygon)+1
1239        assert len(y) == len(polygon)+1
1240        assert x[0] == 1
1241        assert x[1] == 1
1242        assert x[2] == 100
1243        assert x[3] == 1
1244        assert x[4] == 3
1245        assert y[0] == 5
1246        assert y[1] == 1
1247        assert y[2] == 10
1248        assert y[3] == 10
1249        assert y[4] == 6
1250
1251    # Disabled   
1252    def xtest_plot_polygons(self):
1253       
1254        import os
1255       
1256        #Simplest case: Polygon is the unit square
1257        polygon1 = [[0,0], [1,0], [1,1], [0,1]]
1258        polygon2 = [[1,1], [2,1], [3,2], [2,2]]
1259        v = plot_polygons([polygon1, polygon2],'test1')
1260        assert len(v) == 4
1261        assert v[0] == 0
1262        assert v[1] == 3
1263        assert v[2] == 0
1264        assert v[3] == 2
1265
1266        #Another case
1267        polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
1268        v = plot_polygons([polygon2,polygon3],'test2')
1269        assert len(v) == 4
1270        assert v[0] == 1
1271        assert v[1] == 100
1272        assert v[2] == 1
1273        assert v[3] == 10
1274
1275        os.remove('test1.png')
1276        os.remove('test2.png')
1277
1278       
1279    def test_inside_polygon_geospatial(self):
1280
1281
1282        polygon_absolute = [[0,0], [1,0], [1,1], [0,1]]
1283        poly_geo_ref = Geo_reference(57,100,100)
1284       
1285
1286
1287
1288        #Simplest case: Polygon is the unit square
1289        polygon_absolute = [[0,0], [1,0], [1,1], [0,1]]
1290        poly_geo_ref = Geo_reference(57,100,100)
1291        polygon = poly_geo_ref.change_points_geo_ref(polygon_absolute)
1292        poly_spatial = Geospatial_data(polygon,
1293                                       geo_reference=poly_geo_ref)
1294       
1295        points_absolute = (0.5, 0.5)
1296        points_geo_ref = Geo_reference(57,78,-56)
1297        points = points_geo_ref.change_points_geo_ref(points_absolute)
1298        points_spatial = Geospatial_data(points,
1299                                         geo_reference=points_geo_ref) 
1300       
1301        assert is_inside_polygon(points_absolute, polygon_absolute)
1302        assert is_inside_polygon(ensure_numeric(points_absolute),
1303                                 ensure_numeric(polygon_absolute))
1304        assert is_inside_polygon(points_absolute, poly_spatial)
1305        assert is_inside_polygon(points_spatial, poly_spatial)
1306        assert is_inside_polygon(points_spatial, polygon_absolute)
1307
1308        assert is_inside_polygon(points_absolute, polygon_absolute)
1309
1310
1311    def NOtest_decimate_polygon(self):
1312
1313        polygon = [[0,0], [10,10], [15,5], [20, 10],
1314                   [25,0], [30,10], [40,-10], [35, -5]]
1315
1316        #plot_polygons([polygon], figname='test')
1317       
1318        dpoly = decimate_polygon(polygon, factor=2)
1319
1320        print dpoly
1321       
1322        assert len(dpoly)*2==len(polygon)
1323
1324        minx = maxx = polygon[0][0]
1325        miny = maxy = polygon[0][1]
1326        for point in polygon[1:]:
1327            x, y = point
1328           
1329            if x < minx: minx = x
1330            if x > maxx: maxx = x
1331            if y < miny: miny = y
1332            if y > maxy: maxy = y
1333           
1334
1335        assert [minx, miny] in polygon
1336        print minx, maxy
1337        assert [minx, maxy] in polygon
1338        assert [maxx, miny] in polygon
1339        assert [maxx, maxy] in polygon               
1340       
1341
1342       
1343#-------------------------------------------------------------
1344if __name__ == "__main__":
1345    suite = unittest.makeSuite(Test_Polygon,'test')
1346    #suite = unittest.makeSuite(Test_Polygon,'test_inside_polygon_geo_ref')
1347    runner = unittest.TextTestRunner()
1348    runner.run(suite)
1349
1350
1351
1352
Note: See TracBrowser for help on using the repository browser.