Context Navigation

source: trunk/anuga_core/source/anuga/geometry/test_polygon.py @ 8019

Last change on this file since 8019 was 7866, checked in by James Hudson, 15 years ago
More swb tests passing. Cleaned up some pylint errors.
File size: 66.9 KB

Line
1	#!/usr/bin/env python
2
3	""" Test suite to test polygon functionality. """
4
5	import unittest
6	import numpy as num
7	from anuga.utilities.numerical_tools import ensure_numeric
8	from anuga.utilities.system_tools import get_pathname_from_package
9
10	from polygon import _poly_xy, separate_points_by_polygon, \
11	populate_polygon, polygon_area, is_inside_polygon, \
12	read_polygon, point_on_line, point_in_polygon, \
13	plot_polygons, outside_polygon, is_outside_polygon, \
14	intersection, is_complex, is_inside_triangle, \
15	interpolate_polyline, inside_polygon, in_and_outside_polygon
16
17	from polygon_function import Polygon_function
18	from anuga.coordinate_transforms.geo_reference import Geo_reference
19	from anuga.geospatial_data.geospatial_data import Geospatial_data
20
21
22	def test_function(x, y):
23	return x+y
24
25
26	class Test_Polygon(unittest.TestCase):
27	def setUp(self):
28	pass
29
30	def tearDown(self):
31	pass
32
33	def test_that_C_extension_compiles(self):
34	FN = 'polygon_ext.c'
35	try:
36	import polygon_ext
37	except:
38	from compile import compile
39
40	try:
41	compile(FN)
42	except:
43	raise Exception, 'Could not compile %s' %FN
44	else:
45	import polygon_ext
46
47	# Polygon stuff
48	def test_polygon_function_constants(self):
49	p1 = [[0,0], [10,0], [10,10], [0,10]]
50	p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
51
52	f = Polygon_function([(p1, 1.0)])
53	z = f([5, 5, 27, 35], [5, 9, 8, -5]) # Two first inside p1
54	assert num.allclose(z, [1, 1, 0, 0])
55
56	f = Polygon_function([(p2, 2.0)])
57	z = f([5, 5, 27, 35], [5, 9, 8, -5]) # First and last inside p2
58	assert num.allclose(z, [2, 0, 0, 2])
59
60	# Combined
61	f = Polygon_function([(p1, 1.0), (p2, 2.0)])
62	z = f([5, 5, 27, 35], [5, 9, 8, -5])
63	assert num.allclose(z, [2, 1, 0, 2])
64
65	def test_polygon_function_csvfile(self):
66	from os import sep, getenv
67
68	# Get path where this test is run
69	path = get_pathname_from_package('anuga.utilities')
70
71	# Form absolute filename and read
72	filename = path + sep + 'mainland_only.csv'
73	p1 = read_polygon(filename)
74
75	f = Polygon_function([(p1, 10.0)])
76	z = f([430000, 480000], [490000, 7720000]) # 1st outside, 2nd inside
77
78	assert num.allclose(z, [0, 10])
79
80	def test_polygon_function_georef(self):
81	"""Check that georeferencing works"""
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 num.allclose(z, [1, 1, 0, 0])
96
97	f = Polygon_function([(p2, 2.0)], geo_reference=geo)
98	z = f([5, 5, 27, 35], [5, 9, 8, -5]) # First and last inside p2
99	assert num.allclose(z, [2, 0, 0, 2])
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 num.allclose(z, [2, 1, 0, 2])
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 num.allclose(z, [2, 1, 0, 2])
110
111	def test_polygon_function_callable(self):
112	"""Check that values passed into Polygon_function can be callable."""
113
114	p1 = [[0,0], [10,0], [10,10], [0,10]]
115	p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
116
117	f = Polygon_function([(p1, test_function)])
118	z = f([5, 5, 27, 35], [5, 9, 8, -5]) # Two first inside p1
119	assert num.allclose(z, [10, 14, 0, 0])
120
121	# Combined
122	f = Polygon_function([(p1, test_function), (p2, 2.0)])
123	z = f([5, 5, 27, 35], [5, 9, 8, -5])
124	assert num.allclose(z, [2, 14, 0, 2])
125
126	# Combined w default
127	f = Polygon_function([(p1, test_function), (p2, 2.0)], default = 3.14)
128	z = f([5, 5, 27, 35], [5, 9, 8, -5])
129	assert num.allclose(z, [2, 14, 3.14, 2])
130
131	# Combined w default func
132	f = Polygon_function([(p1, test_function), (p2, 2.0)],
133	default=test_function)
134	z = f([5, 5, 27, 35], [5, 9, 8, -5])
135	assert num.allclose(z, [2, 14, 35, 2])
136
137	def test_point_on_line(self):
138	# Endpoints first
139	assert point_on_line([0,0], [[0,0], [1,0]])
140	assert point_on_line([1,0], [[0,0], [1,0]])
141
142	# Endpoints first - non-simple
143	assert point_on_line([-0.1,0.0], [[-0.1,0.0], [1.5,0.6]])
144	assert point_on_line([1.5,0.6], [[-0.1,0.0], [1.5,0.6]])
145
146	# Then points on line
147	assert point_on_line([0.5,0], [[0,0], [1,0]])
148	assert point_on_line([0,0.5], [[0,1], [0,0]])
149	assert point_on_line([1,0.5], [[1,1], [1,0]])
150	assert point_on_line([0.5,0.5], [[0,0], [1,1]])
151
152	# Then points not on line
153	assert not point_on_line([0.5,0], [[0,1], [1,1]])
154	assert not point_on_line([0,0.5], [[0,0], [1,1]])
155
156	# From real example that failed
157	assert not point_on_line([40,50], [[40,20], [40,40]])
158
159	# From real example that failed
160	assert not point_on_line([40,19], [[40,20], [40,40]])
161
162	# Degenerate line
163	assert point_on_line([40,19], [[40,19], [40,19]])
164	assert not point_on_line([40,19], [[40,40], [40,40]])
165
166	def test_is_inside_polygon_main(self):
167	# Simplest case: Polygon is the unit square
168	polygon = [[0,0], [1,0], [1,1], [0,1]]
169
170	assert is_inside_polygon((0.5, 0.5), polygon)
171	assert not is_inside_polygon((0.5, 1.5), polygon)
172	assert not is_inside_polygon((0.5, -0.5), polygon)
173	assert not is_inside_polygon((-0.5, 0.5), polygon)
174	assert not is_inside_polygon((1.5, 0.5), polygon)
175
176	# Try point on borders
177	assert is_inside_polygon((1., 0.5), polygon, closed=True)
178	assert is_inside_polygon((0.5, 1.), polygon, closed=True)
179	assert is_inside_polygon((0., 0.5), polygon, closed=True)
180	assert is_inside_polygon((0.5, 0.), polygon, closed=True)
181
182	assert not is_inside_polygon((0.5, 1.), polygon, closed=False)
183	assert not is_inside_polygon((0., 0.5), polygon, closed=False)
184	assert not is_inside_polygon((0.5, 0.), polygon, closed=False)
185	assert not is_inside_polygon((1., 0.5), polygon, closed=False)
186
187	def test_inside_polygon_main(self):
188	"""test_is_inside_polygon
189
190	Test fast version of of is_inside_polygon
191	"""
192
193	# Simplest case: Polygon is the unit square
194	polygon = [[0,0], [1,0], [1,1], [0,1]]
195
196	assert is_inside_polygon( (0.5, 0.5), polygon )
197	assert not is_inside_polygon( (0.5, 1.5), polygon )
198	assert not is_inside_polygon( (0.5, -0.5), polygon )
199	assert not is_inside_polygon( (-0.5, 0.5), polygon )
200	assert not is_inside_polygon( (1.5, 0.5), polygon )
201
202	# Try point on borders
203	assert is_inside_polygon( (1., 0.5), polygon, closed=True)
204	assert is_inside_polygon( (0.5, 1), polygon, closed=True)
205	assert is_inside_polygon( (0., 0.5), polygon, closed=True)
206	assert is_inside_polygon( (0.5, 0.), polygon, closed=True)
207
208	assert not is_inside_polygon( (0.5, 1), polygon, closed=False)
209	assert not is_inside_polygon( (0., 0.5), polygon, closed=False)
210	assert not is_inside_polygon( (0.5, 0.), polygon, closed=False)
211	assert not is_inside_polygon( (1., 0.5), polygon, closed=False)
212
213
214	def test_inside_polygon_main(self):
215	# Simplest case: Polygon is the unit square
216	polygon = [[0,0], [1,0], [1,1], [0,1]]
217
218	# From real example (that failed)
219	polygon = [[20,20], [40,20], [40,40], [20,40]]
220	points = [[40, 50]]
221	res = inside_polygon(points, polygon)
222	assert len(res) == 0
223
224	polygon = [[20,20], [40,20], [40,40], [20,40]]
225	points = [[25, 25], [30, 20], [40, 50], [90, 20], [40, 90]]
226	res = inside_polygon(points, polygon)
227	assert len(res) == 2
228	assert num.allclose(res, [0,1])
229
230	# More convoluted and non convex polygon
231	polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
232	assert is_inside_polygon( (0.5, 0.5), polygon )
233	assert is_inside_polygon( (1, -0.5), polygon )
234	assert is_inside_polygon( (1.5, 0), polygon )
235
236	assert not is_inside_polygon( (0.5, 1.5), polygon )
237	assert not is_inside_polygon( (0.5, -0.5), polygon )
238
239	assert is_inside_polygon( (0.5, 0.5), polygon )
240	assert is_inside_polygon( (1, -0.5), polygon )
241	assert is_inside_polygon( (1.5, 0), polygon )
242
243	assert not is_inside_polygon( (0.5, 1.5), polygon )
244	assert not is_inside_polygon( (0.5, -0.5), polygon )
245
246	# Very convoluted polygon
247	polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
248	assert is_inside_polygon((5, 5), polygon)
249	assert is_inside_polygon((17, 7), polygon)
250	assert is_inside_polygon((27, 2), polygon)
251	assert is_inside_polygon((35, -5), polygon)
252	assert not is_inside_polygon((15, 7), polygon)
253	assert not is_inside_polygon((24, 3), polygon)
254	assert not is_inside_polygon((25, -10), polygon)
255
256	# Another combination (that failed)
257	polygon = [[0,0], [10,0], [10,10], [0,10]]
258	assert is_inside_polygon((5, 5), polygon)
259	assert is_inside_polygon((7, 7), polygon)
260	assert not is_inside_polygon((-17, 7), polygon)
261	assert not is_inside_polygon((7, 17), polygon)
262	assert not is_inside_polygon((17, 7), polygon)
263	assert not is_inside_polygon((27, 8), polygon)
264	assert not is_inside_polygon((35, -5), polygon)
265
266	# Multiple polygons
267	polygon = [[0,0], [1,0], [1,1], [0,1], [0,0], [10,10],
268	[11,10], [11,11], [10,11], [10,10]]
269	assert is_inside_polygon((0.5, 0.5), polygon)
270	assert is_inside_polygon((10.5, 10.5), polygon)
271
272	# FIXME: Fails if point is 5.5, 5.5
273	assert not is_inside_polygon((0, 5.5), polygon)
274
275	# Polygon with a hole
276	polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
277	[0,0], [1,0], [1,1], [0,1], [0,0]]
278
279	assert is_inside_polygon((0, -0.5), polygon)
280	assert not is_inside_polygon((0.5, 0.5), polygon)
281
282	def test_duplicate_points_being_ok(self):
283	# Simplest case: Polygon is the unit square
284	polygon = [[0,0], [1,0], [1,0], [1,0], [1,1], [0,1], [0,0]]
285
286	assert is_inside_polygon((0.5, 0.5), polygon)
287	assert not is_inside_polygon((0.5, 1.5), polygon)
288	assert not is_inside_polygon((0.5, -0.5), polygon)
289	assert not is_inside_polygon((-0.5, 0.5), polygon)
290	assert not is_inside_polygon((1.5, 0.5), polygon)
291
292	# Try point on borders
293	assert is_inside_polygon((1., 0.5), polygon, closed=True)
294	assert is_inside_polygon((0.5, 1), polygon, closed=True)
295	assert is_inside_polygon((0., 0.5), polygon, closed=True)
296	assert is_inside_polygon((0.5, 0.), polygon, closed=True)
297
298	assert not is_inside_polygon((0.5, 1), polygon, closed=False)
299	assert not is_inside_polygon((0., 0.5), polygon, closed=False)
300	assert not is_inside_polygon((0.5, 0.), polygon, closed=False)
301	assert not is_inside_polygon((1., 0.5), polygon, closed=False)
302
303	# From real example
304	polygon = [[20,20], [40,20], [40,40], [20,40]]
305	points = [[40, 50]]
306	res = inside_polygon(points, polygon)
307	assert len(res) == 0
308
309	def test_inside_polygon_vector_version(self):
310	# Now try the vector formulation returning indices
311	polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
312	points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
313	res = inside_polygon( points, polygon, verbose=False )
314	assert num.allclose(res, [0,1,2])
315
316	def test_outside_polygon(self):
317	# unit square
318	U = [[0,0], [1,0], [1,1], [0,1]]
319
320	# evaluate to False as the point 0.5, 0.5 is inside the unit square
321	assert not is_outside_polygon([0.5, 0.5], U)
322
323	# evaluate to True as the point 1.5, 0.5 is outside the unit square
324	assert is_outside_polygon([1.5, 0.5], U)
325
326	indices = outside_polygon([[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
327	assert num.allclose(indices, [1])
328
329	# One more test of vector formulation returning indices
330	polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
331	points = [[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	assert num.allclose(res, [3, 4])
334
335	polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
336	points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
337	[0.5, 1.5], [0.5, -0.5]]
338	res = outside_polygon(points, polygon)
339
340	assert num.allclose(res, [0, 4, 5])
341
342	def test_outside_polygon2(self):
343	# unit square
344	U = [[0,0], [1,0], [1,1], [0,1]]
345
346	# evaluate to False as the point 0.5, 1.0 is inside the unit square
347	assert not outside_polygon([0.5, 1.0], U, closed = True)
348
349	# evaluate to True as the point 0.5, 1.0 is outside the unit square
350	assert is_outside_polygon([0.5, 1.0], U, closed = False)
351
352	def test_all_outside_polygon(self):
353	"""Test case where all points are outside poly"""
354
355	# unit square
356	U = [[0,0], [1,0], [1,1], [0,1]]
357
358	points = [[2,2], [1,3], [-1,1], [0,2]] # All outside
359
360	indices, count = separate_points_by_polygon(points, U)
361	assert count == 0 # None inside
362	assert num.allclose(indices, [3, 2, 1, 0])
363
364	indices = outside_polygon(points, U, closed = True)
365	assert num.allclose(indices, [0, 1, 2, 3])
366
367	indices = inside_polygon(points, U, closed = True)
368	assert num.allclose(indices, [])
369
370	def test_all_inside_polygon(self):
371	"""Test case where all points are inside poly"""
372
373	# unit square
374	U = [[0,0], [1,0], [1,1], [0,1]]
375
376	points = [[0.5,0.5], [0.2,0.3], [0,0.5]] # All inside (or on edge)
377
378	indices, count = separate_points_by_polygon(points, U)
379	assert count == 3 # All inside
380	assert num.allclose(indices, [0, 1, 2])
381
382	indices = outside_polygon(points, U, closed=True)
383	assert num.allclose(indices, [])
384
385	indices = inside_polygon(points, U, closed=True)
386	assert num.allclose(indices, [0, 1, 2])
387
388
389	def test_separate_points_by_polygon(self):
390	# unit square
391	U = [[0,0], [1,0], [1,1], [0,1]]
392
393	indices, count = separate_points_by_polygon([[0.5, 0.5],
394	[1, -0.5],
395	[0.3, 0.2]],
396	U)
397	assert num.allclose( indices, [0, 2, 1] )
398	assert count == 2
399
400	# One more test of vector formulation returning indices
401	polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
402	points = [[0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
403	res, count = separate_points_by_polygon(points, polygon)
404
405	assert num.allclose(res, [0, 1, 2, 4, 3])
406	assert count == 3
407
408	polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
409	points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], \
410	[0.5, -0.5]]
411	res, count = separate_points_by_polygon( points, polygon )
412
413	assert num.allclose( res, [1,2,3,5,4,0] )
414	assert count == 3
415
416
417	def test_is_inside_triangle(self):
418	# Simplest case:
419	triangle = [[0, 0], [1, 0], [0.5, 1]]
420
421	assert is_inside_triangle((0.5, 0.5), triangle)
422	assert is_inside_triangle((0.9, 0.1), triangle)
423	assert not is_inside_triangle((0.5, 1.5), triangle)
424	assert not is_inside_triangle((0.5, -0.5), triangle)
425	assert not is_inside_triangle((-0.5, 0.5), triangle)
426	assert not is_inside_triangle((1.5, 0.5), triangle)
427
428	# Try point on borders
429	assert is_inside_triangle((0.5, 0), triangle, closed=True)
430	assert is_inside_triangle((1, 0), triangle, closed=True)
431
432	assert not is_inside_triangle((0.5, 0), triangle, closed=False)
433	assert not is_inside_triangle((1, 0), triangle, closed=False)
434
435	# Try vertices
436	for P in triangle:
437	assert is_inside_triangle(P, triangle, closed=True)
438	assert not is_inside_triangle(P, triangle, closed=False)
439
440	# Slightly different
441	triangle = [[0, 0.1], [1, -0.2], [0.5, 1]]
442	assert is_inside_triangle((0.5, 0.5), triangle)
443	assert is_inside_triangle((0.4, 0.1), triangle)
444	assert not is_inside_triangle((1, 1), triangle)
445
446	# Try vertices
447	for P in triangle:
448	assert is_inside_triangle(P, triangle, closed=True)
449	assert not is_inside_triangle(P, triangle, closed=False)
450
451
452	def test_populate_polygon(self):
453	polygon = [[0,0], [1,0], [1,1], [0,1]]
454	points = populate_polygon(polygon, 5)
455
456	assert len(points) == 5
457	for point in points:
458	assert is_inside_polygon(point, polygon)
459
460	assert is_inside_polygon(point, polygon)
461
462
463	# Very convoluted polygon
464	polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
465	points = populate_polygon(polygon, 5)
466	assert len(points) == 5
467	for point in points:
468	assert is_inside_polygon(point, polygon)
469	assert is_inside_polygon(point, polygon)
470
471
472	def test_populate_polygon_with_exclude(self):
473	polygon = [[0,0], [1,0], [1,1], [0,1]]
474	ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] # SW quarter
475	points = populate_polygon(polygon, 5, exclude=[ex_poly])
476
477	assert len(points) == 5
478	for point in points:
479	assert is_inside_polygon(point, polygon)
480	assert not is_inside_polygon(point, ex_poly)
481
482	# overlap
483	polygon = [[0,0], [1,0], [1,1], [0,1]]
484	ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
485	points = populate_polygon(polygon, 5, exclude=[ex_poly])
486
487	assert len(points) == 5
488	for point in points:
489	assert is_inside_polygon(point, polygon)
490	assert not is_inside_polygon(point, ex_poly)
491
492	# Multiple
493	polygon = [[0,0], [1,0], [1,1], [0,1]]
494	ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] # SW quarter
495	ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] # NE quarter
496
497	points = populate_polygon(polygon, 20, exclude=[ex_poly1, ex_poly2])
498
499	assert len(points) == 20
500	for point in points:
501	assert is_inside_polygon(point, polygon)
502	assert not is_inside_polygon(point, ex_poly1)
503	assert not is_inside_polygon(point, ex_poly2)
504
505	# Very convoluted polygon
506	polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
507	ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
508	points = populate_polygon(polygon, 20, exclude=[ex_poly])
509
510	assert len(points) == 20
511	for point in points:
512	assert is_inside_polygon(point, polygon)
513	assert not is_inside_polygon(point, ex_poly), '%s' % str(point)
514
515	def test_populate_polygon_with_exclude2(self):
516	min_outer = 0
517	max_outer = 1000
518	polygon_outer = [[min_outer,min_outer], [max_outer,min_outer],
519	[max_outer,max_outer], [min_outer,max_outer]]
520
521	delta = 10
522	min_inner1 = min_outer + delta
523	max_inner1 = max_outer - delta
524	inner1_polygon = [[min_inner1,min_inner1], [max_inner1,min_inner1],
525	[max_inner1,max_inner1], [min_inner1,max_inner1]]
526
527	density_inner2 = 1000
528	min_inner2 = min_outer + 2*delta
529	max_inner2 = max_outer - 2*delta
530	inner2_polygon = [[min_inner2,min_inner2], [max_inner2,min_inner2],
531	[max_inner2,max_inner2], [min_inner2,max_inner2]]
532
533	points = populate_polygon(polygon_outer, 20,
534	exclude=[inner1_polygon, inner2_polygon])
535
536	assert len(points) == 20
537	for point in points:
538	assert is_inside_polygon(point, polygon_outer)
539	assert not is_inside_polygon(point, inner1_polygon)
540	assert not is_inside_polygon(point, inner2_polygon)
541
542	# Very convoluted polygon
543	polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
544	ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
545	points = populate_polygon(polygon, 20, exclude=[ex_poly])
546
547	assert len(points) == 20
548	for point in points:
549	assert is_inside_polygon(point, polygon)
550	assert not is_inside_polygon(point, ex_poly), '%s' % str(point)
551
552	def test_point_in_polygon(self):
553	polygon = [[0,0], [1,0], [1,1], [0,1]]
554	point = point_in_polygon(polygon)
555	assert is_inside_polygon(point, polygon)
556
557	# this may get into a vicious loop
558	# polygon = [[1e32,1e54], [1,0], [1,1], [0,1]]
559
560	polygon = [[1e15,1e7], [1,0], [1,1], [0,1]]
561	point = point_in_polygon(polygon)
562	assert is_inside_polygon(point, polygon)
563
564	polygon = [[0,0], [1,0], [1,1], [1e8,1e8]]
565	point = point_in_polygon(polygon)
566	assert is_inside_polygon(point, polygon)
567
568	polygon = [[1e32,1e54], [-1e32,1e54], [1e32,-1e54]]
569	point = point_in_polygon(polygon)
570	assert is_inside_polygon(point, polygon)
571
572	polygon = [[1e18,1e15], [1,0], [0,1]]
573	point = point_in_polygon(polygon)
574	assert is_inside_polygon(point, polygon)
575
576	def test_in_and_outside_polygon_main(self):
577	# Simplest case: Polygon is the unit square
578	polygon = [[0,0], [1,0], [1,1], [0,1]]
579
580	inside, outside = in_and_outside_polygon((0.5, 0.5), polygon)
581	assert inside[0] == 0
582	assert len(outside) == 0
583
584	inside, outside = in_and_outside_polygon((1., 0.5), polygon,
585	closed=True)
586	assert inside[0] == 0
587	assert len(outside) == 0
588
589	inside, outside = in_and_outside_polygon((1., 0.5), polygon,
590	closed=False)
591	assert len(inside) == 0
592	assert outside[0] == 0
593
594	points = [(1., 0.25),(1., 0.75)]
595	inside, outside = in_and_outside_polygon(points, polygon, closed=True)
596	assert num.alltrue(inside == [0,1])
597	assert len(outside) == 0
598
599	inside, outside = in_and_outside_polygon(points, polygon, closed=False)
600	assert len(inside) == 0
601	assert num.alltrue(outside == [0,1])
602
603	points = [(100., 0.25), (0.5, 0.5) ]
604	inside, outside = in_and_outside_polygon(points, polygon)
605	assert num.alltrue(inside == [1])
606	assert outside[0] == 0
607
608	points = [(100., 0.25),(0.5, 0.5), (39,20), (0.6,0.7),(56,43),(67,90)]
609	inside, outside = in_and_outside_polygon(points, polygon)
610	assert num.alltrue(inside == [1, 3])
611	assert num.alltrue(outside == [0, 2, 4, 5])
612
613	def test_intersection1(self):
614	line0 = [[-1,0], [1,0]]
615	line1 = [[0,-1], [0,1]]
616
617	status, value = intersection(line0, line1)
618	assert status == 1
619	assert num.allclose(value, [0.0, 0.0])
620
621	def test_intersection2(self):
622	line0 = [[0,0], [24,12]]
623	line1 = [[0,12], [24,0]]
624
625	status, value = intersection(line0, line1)
626	assert status == 1
627	assert num.allclose(value, [12.0, 6.0])
628
629	# Swap direction of one line
630	line1 = [[24,0], [0,12]]
631
632	status, value = intersection(line0, line1)
633	assert status == 1
634	assert num.allclose(value, [12.0, 6.0])
635
636	# Swap order of lines
637	status, value = intersection(line1, line0)
638	assert status == 1
639	assert num.allclose(value, [12.0, 6.0])
640
641	def test_intersection3(self):
642	line0 = [[0,0], [24,12]]
643	line1 = [[0,17], [24,0]]
644
645	status, value = intersection(line0, line1)
646	assert status == 1
647	assert num.allclose(value, [14.068965517, 7.0344827586])
648
649	# Swap direction of one line
650	line1 = [[24,0], [0,17]]
651
652	status, value = intersection(line0, line1)
653	assert status == 1
654	assert num.allclose(value, [14.068965517, 7.0344827586])
655
656	# Swap order of lines
657	status, value = intersection(line1, line0)
658	assert status == 1
659	assert num.allclose(value, [14.068965517, 7.0344827586])
660
661	def test_intersection_endpoints(self):
662	"""test_intersection_endpoints(self):
663
664	Test that coinciding endpoints are picked up
665	"""
666
667	line0 = [[0,0], [1,1]]
668	line1 = [[1,1], [2,1]]
669
670	status, value = intersection(line0, line1)
671	assert status == 1
672	assert num.allclose(value, [1.0, 1.0])
673
674	line0 = [[1,1], [2,0]]
675	line1 = [[1,1], [2,1]]
676
677	status, value = intersection(line0, line1)
678	assert status == 1
679	assert num.allclose(value, [1.0, 1.0])
680
681	# This function is a helper function for
682	# the test_intersection_bug_20081110_?() set of tests.
683	# This function tests all parallel line cases for 4 collinear points.
684	# This function should never be run directly by the unittest code.
685	def helper_test_parallel_intersection_code(self, P1, P2, P3, P4):
686	# lines in same direction, no overlap
687	# 0: ---->----
688	# 1: --------->-----------
689	line0 = [P1, P2]
690	line1 = [P3, P4]
691	status, value = intersection(line0, line1)
692	self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
693	(str(status), str(value)))
694	self.failUnless(value is None, 'Expected value of None, got %s' %
695	str(value))
696
697	# lines in same direction, no overlap
698	# 0: ----<----
699	# 1: ---------<-----------
700	line0 = [P2, P1]
701	line1 = [P4, P3]
702	status, value = intersection(line0, line1)
703	self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
704	(str(status), str(value)))
705	self.failUnless(value is None, 'Expected value of None, got %s' %
706	str(value))
707
708	# lines in opposite direction, no overlap
709	# 0: ----<----
710	# 1: --------->-----------
711	line0 = [P2, P1]
712	line1 = [P3, P4]
713	status, value = intersection(line0, line1)
714	self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
715	(str(status), str(value)))
716	self.failUnless(value is None, 'Expected value of None, got %s' %
717	str(value))
718
719	# lines in opposite direction, no overlap
720	# 0: ---->----
721	# 1: ---------<-----------
722	line0 = [P1, P2]
723	line1 = [P4, P3]
724	status, value = intersection(line0, line1)
725	self.failIf(status!=3, 'Expected status 3, got status=%s, value=%s' %
726	(str(status), str(value)))
727	self.failUnless(value is None, 'Expected value of None, got %s' %
728	str(value))
729
730	# ----------------------------------------------------------------------
731
732	# line0 fully within line1, same direction
733	# 0: ---->----
734	# 1: --------->-----------
735	# value should be line0:
736	# ---->----
737	line0 = [P2, P3]
738	line1 = [P1, P4]
739	status, value = intersection(line0, line1)
740	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
741	(str(status), str(value)))
742	self.failUnless(num.allclose(value, line0))
743
744	# line0 fully within line1, same direction
745	# 0: ----<----
746	# 1: ---------<-----------
747	# value should be line0:
748	# ----<----
749	line0 = [P3, P2]
750	line1 = [P4, P1]
751	status, value = intersection(line0, line1)
752	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
753	(str(status), str(value)))
754	self.failUnless(num.allclose(value, line0))
755
756	# line0 fully within line1, opposite direction
757	# 0: ----<----
758	# 1: --------->-----------
759	# value should be line0:
760	# ----<----
761	line0 = [P3, P2]
762	line1 = [P1, P4]
763	status, value = intersection(line0, line1)
764	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
765	(str(status), str(value)))
766	self.failUnless(num.allclose(value, line0))
767
768	# line0 fully within line1, opposite direction
769	# 0: ---->----
770	# 1: ---------<-----------
771	# value should be line0:
772	# ---->----
773	line0 = [P2, P3]
774	line1 = [P4, P1]
775	status, value = intersection(line0, line1)
776	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
777	(str(status), str(value)))
778	self.failUnless(num.allclose(value, line0))
779
780	# ----------------------------------------------------------------------
781
782	# line1 fully within line0, same direction
783	# 0: --------->-----------
784	# 1: ---->----
785	# value should be line1:
786	# ---->----
787	line0 = [P1, P4]
788	line1 = [P2, P3]
789	status, value = intersection(line0, line1)
790	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
791	(str(status), str(value)))
792	self.failUnless(num.allclose(value, line1))
793
794	# line1 fully within line0, same direction
795	# 0: ---------<-----------
796	# 1: ----<----
797	# value should be line1:
798	# ----<----
799	line0 = [P4, P1]
800	line1 = [P3, P2]
801	status, value = intersection(line0, line1)
802	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
803	(str(status), str(value)))
804	self.failUnless(num.allclose(value, line1))
805
806	# line1 fully within line0, opposite direction
807	# 0: ---------<-----------
808	# 1: ---->----
809	# value should be line1:
810	# ---->----
811	line0 = [P4, P1]
812	line1 = [P2, P3]
813	status, value = intersection(line0, line1)
814	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
815	(str(status), str(value)))
816	self.failUnless(num.allclose(value, line1))
817
818	# line1 fully within line0, opposite direction
819	# 0: --------->-----------
820	# 1: ----<----
821	# value should be line1:
822	# ----<----
823	line0 = [P1, P4]
824	line1 = [P3, P2]
825	status, value = intersection(line0, line1)
826	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
827	(str(status), str(value)))
828	self.failUnless(num.allclose(value, line1))
829
830	# ----------------------------------------------------------------------
831
832	# line in same direction, partial overlap
833	# 0: ----->-----
834	# 1: ------->--------
835	# value should be segment line1_start->line0_end:
836	# --->----
837	line0 = [P1, P3]
838	line1 = [P2, P4]
839	status, value = intersection(line0, line1)
840	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
841	(str(status), str(value)))
842	self.failUnless(num.allclose(value, [line1[0],line0[1]]))
843
844	# line in same direction, partial overlap
845	# 0: -----<-----
846	# 1: -------<--------
847	# value should be segment line0_start->line1_end:
848	# ---<----
849	line0 = [P3, P1]
850	line1 = [P4, P2]
851	status, value = intersection(line0, line1)
852	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
853	(str(status), str(value)))
854	self.failUnless(num.allclose(value, [line0[0],line1[1]]))
855
856	# line in opposite direction, partial overlap
857	# 0: -----<-----
858	# 1: ------->--------
859	# value should be segment line0_start->line1_start:
860	# ---<----
861	line0 = [P3, P1]
862	line1 = [P2, P4]
863	status, value = intersection(line0, line1)
864	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
865	(str(status), str(value)))
866	self.failUnless(num.allclose(value, [line0[0],line1[0]]))
867
868	# line in opposite direction, partial overlap
869	# 0: ----->-----
870	# 1: -------<--------
871	# value should be segment line1_end->line0_end:
872	# --->----
873	line0 = [P1, P3]
874	line1 = [P4, P2]
875	status, value = intersection(line0, line1)
876	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
877	(str(status), str(value)))
878	self.failUnless(num.allclose(value, [line1[1],line0[1]]))
879
880	# ----------------------------------------------------------------------
881
882	# line in same direction, partial overlap
883	# 0: ------>------
884	# 1: ------>------
885	# value should be segment line0_start->line1_end:
886	# --->----
887	line0 = [P2, P4]
888	line1 = [P1, P3]
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(num.allclose(value, [line0[0],line1[1]]))
893
894	# line in same direction, partial overlap
895	# 0: ------<------
896	# 1: ------<------
897	# value should be segment line1_start->line0_end:
898	# ----<-----
899	line0 = [P4, P2]
900	line1 = [P3, P1]
901	status, value = intersection(line0, line1)
902	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
903	(str(status), str(value)))
904	self.failUnless(num.allclose(value, [line1[0],line0[1]]))
905
906	# line in opposite direction, partial overlap
907	# 0: ------<------
908	# 1: ----->------
909	# value should be segment line1_end->line0_end:
910	# --->----
911	line0 = [P4, P2]
912	line1 = [P1, P3]
913	status, value = intersection(line0, line1)
914	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
915	(str(status), str(value)))
916	self.failUnless(num.allclose(value, [line1[1],line0[1]]))
917
918	# line in opposite direction, partial overlap
919	# 0: ------>------
920	# 1: -----<------
921	# value should be segment line0_start->line1_start:
922	# ---<----
923	line0 = [P2, P4]
924	line1 = [P3, P1]
925	status, value = intersection(line0, line1)
926	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
927	(str(status), str(value)))
928	self.failUnless(num.allclose(value, [line0[0],line1[0]]))
929
930	# ----------------------------------------------------------------------
931
932	# line in same direction, same left point, line1 longer
933	# 0: ----->------
934	# 1: ------->--------
935	# value should be line0:
936	# ----->------
937	line0 = [P1, P3]
938	line1 = [P1, P4]
939	status, value = intersection(line0, line1)
940	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
941	(str(status), str(value)))
942	self.failUnless(num.allclose(value, line0))
943
944	# line in same direction, same left point, line1 longer
945	# 0: -----<------
946	# 1: -------<--------
947	# value should be line0:
948	# -----<------
949	line0 = [P3, P1]
950	line1 = [P4, P1]
951	status, value = intersection(line0, line1)
952	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
953	(str(status), str(value)))
954	self.failUnless(num.allclose(value, line0))
955
956	# line in opposite direction, same left point, line1 longer
957	# 0: ----->------
958	# 1: -------<--------
959	# value should be line0:
960	# ----->------
961	line0 = [P1, P3]
962	line1 = [P4, P1]
963	status, value = intersection(line0, line1)
964	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
965	(str(status), str(value)))
966	self.failUnless(num.allclose(value, line0))
967
968	# line in opposite direction, same start point, line1 longer
969	# 0: -----<------
970	# 1: ------->--------
971	# value should be line0:
972	# -----<------
973	line0 = [P3, P1]
974	line1 = [P1, P4]
975	status, value = intersection(line0, line1)
976	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
977	(str(status), str(value)))
978	self.failUnless(num.allclose(value, line0))
979
980	# ----------------------------------------------------------------------
981
982	# line in same direction, same left point, same right point
983	# 0: ------->--------
984	# 1: ------->--------
985	# value should be line0 or line1:
986	# ------->--------
987	line0 = [P1, P3]
988	line1 = [P1, P3]
989	status, value = intersection(line0, line1)
990	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
991	(str(status), str(value)))
992	self.failUnless(num.allclose(value, line0))
993
994	# line in same direction, same left point, same right point
995	# 0: -------<--------
996	# 1: -------<--------
997	# value should be line0 (or line1):
998	# -------<--------
999	line0 = [P3, P1]
1000	line1 = [P3, P1]
1001	status, value = intersection(line0, line1)
1002	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1003	(str(status), str(value)))
1004	self.failUnless(num.allclose(value, line0))
1005
1006	# line in opposite direction, same left point, same right point
1007	# 0: ------->--------
1008	# 1: -------<--------
1009	# value should be line0:
1010	# ------->--------
1011	line0 = [P1, P3]
1012	line1 = [P3, P1]
1013	status, value = intersection(line0, line1)
1014	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1015	(str(status), str(value)))
1016	self.failUnless(num.allclose(value, line0))
1017
1018	# line in opposite direction, same left point, same right point
1019	# 0: -------<--------
1020	# 1: ------->--------
1021	# value should be line0:
1022	# -------<--------
1023	line0 = [P3, P1]
1024	line1 = [P1, P3]
1025	status, value = intersection(line0, line1)
1026	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1027	(str(status), str(value)))
1028	self.failUnless(num.allclose(value, line0))
1029
1030	# ----------------------------------------------------------------------
1031
1032	# line in same direction, same right point, line1 longer
1033	# 0: ----->------
1034	# 1: ------->--------
1035	# value should be line0:
1036	# ----->------
1037	line0 = [P2, P4]
1038	line1 = [P1, P4]
1039	status, value = intersection(line0, line1)
1040	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1041	(str(status), str(value)))
1042	self.failUnless(num.allclose(value, line0))
1043
1044	# line in same direction, same right point, line1 longer
1045	# 0: -----<------
1046	# 1: -------<--------
1047	# value should be line0:
1048	# -----<------
1049	line0 = [P4, P2]
1050	line1 = [P4, P1]
1051	status, value = intersection(line0, line1)
1052	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1053	(str(status), str(value)))
1054	self.failUnless(num.allclose(value, line0))
1055
1056	# line in opposite direction, same right point, line1 longer
1057	# 0: ----->------
1058	# 1: -------<--------
1059	# value should be line0:
1060	# ----->------
1061	line0 = [P2, P4]
1062	line1 = [P4, P1]
1063	status, value = intersection(line0, line1)
1064	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1065	(str(status), str(value)))
1066	self.failUnless(num.allclose(value, line0))
1067
1068	# line in opposite direction, same right point, line1 longer
1069	# 0: -----<------
1070	# 1: ------->--------
1071	# value should be line0:
1072	# -----<------
1073	line0 = [P4, P2]
1074	line1 = [P1, P4]
1075	status, value = intersection(line0, line1)
1076	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1077	(str(status), str(value)))
1078	self.failUnless(num.allclose(value, line0))
1079
1080	# ----------------------------------------------------------------------
1081
1082	# line in same direction, same left point, line0 longer
1083	# 0: ------->--------
1084	# 1: ----->------
1085	# value should be line1:
1086	# ----->------
1087	line0 = [P1, P4]
1088	line1 = [P1, P3]
1089	status, value = intersection(line0, line1)
1090	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1091	(str(status), str(value)))
1092	self.failUnless(num.allclose(value, line1))
1093
1094	# line in same direction, same left point, line0 longer
1095	# 0: -------<--------
1096	# 1: -----<------
1097	# value should be line1:
1098	# -----<------
1099	line0 = [P4, P1]
1100	line1 = [P3, P1]
1101	status, value = intersection(line0, line1)
1102	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1103	(str(status), str(value)))
1104	self.failUnless(num.allclose(value, line1))
1105
1106	# line in opposite direction, same left point, line0 longer
1107	# 0: ------->--------
1108	# 1: -----<------
1109	# value should be line1:
1110	# -----<------
1111	line0 = [P1, P4]
1112	line1 = [P3, P1]
1113	status, value = intersection(line0, line1)
1114	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1115	(str(status), str(value)))
1116	self.failUnless(num.allclose(value, line1))
1117
1118	# line in opposite direction, same left point, line0 longer
1119	# 0: -------<--------
1120	# 1: ----->------
1121	# value should be line1:
1122	# ----->------
1123	line0 = [P4, P1]
1124	line1 = [P1, P3]
1125	status, value = intersection(line0, line1)
1126	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1127	(str(status), str(value)))
1128	self.failUnless(num.allclose(value, line1))
1129
1130	# ----------------------------------------------------------------------
1131
1132	# line in same direction, same right point, line0 longer
1133	# 0: ------->--------
1134	# 1: ----->------
1135	# value should be line1:
1136	# ----->------
1137	line0 = [P1, P4]
1138	line1 = [P2, P4]
1139	status, value = intersection(line0, line1)
1140	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1141	(str(status), str(value)))
1142	self.failUnless(num.allclose(value, line1))
1143
1144	# line in same direction, same right point, line0 longer
1145	# 0: -------<--------
1146	# 1: -----<------
1147	# value should be line1:
1148	# -----<------
1149	line0 = [P4, P1]
1150	line1 = [P4, P2]
1151	status, value = intersection(line0, line1)
1152	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1153	(str(status), str(value)))
1154	self.failUnless(num.allclose(value, line1))
1155
1156	# line in opposite direction, same right point, line0 longer
1157	# 0: ------->--------
1158	# 1: -----<------
1159	# value should be line1:
1160	# -----<------
1161	line0 = [P1, P4]
1162	line1 = [P4, P2]
1163	status, value = intersection(line0, line1)
1164	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1165	(str(status), str(value)))
1166	self.failUnless(num.allclose(value, line1))
1167
1168	# line in opposite direction, same right point, line0 longer
1169	# 0: -------<--------
1170	# 1: ----->------
1171	# value should be line1:
1172	# ----->------
1173	line0 = [P4, P1]
1174	line1 = [P2, P4]
1175	status, value = intersection(line0, line1)
1176	self.failIf(status!=2, 'Expected status 2, got status=%s, value=%s' %
1177	(str(status), str(value)))
1178	self.failUnless(num.allclose(value, line1))
1179
1180
1181	def test_intersection_bug_20081110_TR(self):
1182	"""test_intersection_bug_20081110(self)
1183
1184	Test all cases in top-right quadrant
1185	"""
1186
1187	# define 4 collinear points in top-right quadrant
1188	# P1---P2---P3---P4
1189	P1 = [1.0, 1.0]
1190	P2 = [2.0, 2.0]
1191	P3 = [3.0, 3.0]
1192	P4 = [4.0, 4.0]
1193
1194	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1195	P1 = [1.0, 1.0+1.0e-9]
1196	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1197	P1 = [1.0, 1.0]
1198	P2 = [2.0, 2.0+1.0e-9]
1199	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1200	P2 = [2.0, 2.0]
1201	P3 = [3.0, 3.0+1.0e-9]
1202	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1203	P3 = [3.0, 3.0]
1204	P4 = [4.0, 4.0+1.0e-9]
1205	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1206
1207	def test_intersection_bug_20081110_TL(self):
1208	"""test_intersection_bug_20081110(self)
1209
1210	Test all cases in top-left quadrant
1211	"""
1212
1213	# define 4 collinear points in top-left quadrant
1214	# P1---P2---P3---P4
1215	P1 = [-1.0, 1.0]
1216	P2 = [-2.0, 2.0]
1217	P3 = [-3.0, 3.0]
1218	P4 = [-4.0, 4.0]
1219
1220	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1221	P1 = [-1.0, 1.0+1.0e-9]
1222	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1223	P1 = [-1.0, 1.0]
1224	P2 = [-2.0, 2.0+1.0e-9]
1225	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1226	P2 = [-2.0, 2.0]
1227	P3 = [-3.0, 3.0+1.0e-9]
1228	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1229	P3 = [-3.0, 3.0]
1230	P4 = [-4.0, 4.0+1.0e-9]
1231	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1232
1233	def test_intersection_bug_20081110_BL(self):
1234	"""test_intersection_bug_20081110(self)
1235
1236	Test all cases in bottom-left quadrant
1237	"""
1238
1239	# define 4 collinear points in bottom-left quadrant
1240	# P1---P2---P3---P4
1241	P1 = [-1.0, -1.0]
1242	P2 = [-2.0, -2.0]
1243	P3 = [-3.0, -3.0]
1244	P4 = [-4.0, -4.0]
1245
1246	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1247	P1 = [-1.0, -1.0+1.0e-9]
1248	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1249	P1 = [-1.0, -1.0]
1250	P2 = [-2.0, -2.0+1.0e-9]
1251	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1252	P2 = [-2.0, -2.0]
1253	P3 = [-3.0, -3.0+1.0e-9]
1254	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1255	P3 = [-3.0, -3.0]
1256	P4 = [-4.0, -4.0+1.0e-9]
1257	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1258
1259	def test_intersection_bug_20081110_BR(self):
1260	"""test_intersection_bug_20081110(self)
1261
1262	Test all cases in bottom-right quadrant
1263	"""
1264
1265	# define 4 collinear points in bottom-right quadrant
1266	# P1---P2---P3---P4
1267	P1 = [1.0, -1.0]
1268	P2 = [2.0, -2.0]
1269	P3 = [3.0, -3.0]
1270	P4 = [4.0, -4.0]
1271
1272	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1273	P1 = [1.0, -1.0+1.0e-9]
1274	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1275	P1 = [1.0, -1.0]
1276	P2 = [2.0, -2.0+1.0e-9]
1277	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1278	P2 = [2.0, -2.0]
1279	P3 = [3.0, -3.0+1.0e-9]
1280	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1281	P3 = [3.0, -3.0]
1282	P4 = [4.0, -4.0+1.0e-9]
1283	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1284
1285	def test_intersection_bug_20081110_TR_TL(self):
1286	"""test_intersection_bug_20081110(self)
1287
1288	Test all cases in top-right & top-left quadrant
1289	"""
1290
1291	# define 4 collinear points, 1 in TL, 3 in TR
1292	# P1-+-P2---P3---P4
1293	P1 = [-3.0, 1.0]
1294	P2 = [ 1.0, 5.0]
1295	P3 = [ 2.0, 6.0]
1296	P4 = [ 3.0, 7.0]
1297	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1298
1299	# define 4 collinear points, 2 in TL, 2 in TR
1300	# P1---P2-+-P3---P4
1301	P1 = [-3.0, 1.0]
1302	P2 = [-2.0, 2.0]
1303	P3 = [ 2.0, 6.0]
1304	P4 = [ 3.0, 7.0]
1305	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1306
1307	# define 4 collinear points, 3 in TL, 1 in TR
1308	# P1---P2---P3-+-P4
1309	P1 = [-3.0, 1.0]
1310	P2 = [-2.0, 2.0]
1311	P3 = [-1.0, 3.0]
1312	P4 = [ 3.0, 7.0]
1313	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1314
1315	def test_intersection_bug_20081110_TR_BL(self):
1316	"""test_intersection_bug_20081110(self)
1317
1318	Test all cases in top-right & bottom-left quadrant
1319	"""
1320
1321	# define 4 collinear points, 1 in BL, 3 in TR
1322	# P1-+-P2---P3---P4
1323	P1 = [-4.0, -3.0]
1324	P2 = [ 1.0, 2.0]
1325	P3 = [ 2.0, 3.0]
1326	P4 = [ 3.0, 4.0]
1327	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1328
1329	# define 4 collinear points, 2 in TL, 2 in TR
1330	# P1---P2-+-P3---P4
1331	P1 = [-4.0, -3.0]
1332	P2 = [-3.0, -2.0]
1333	P3 = [ 2.0, 3.0]
1334	P4 = [ 3.0, 4.0]
1335	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1336
1337	# define 4 collinear points, 3 in TL, 1 in TR
1338	# P1---P2---P3-+-P4
1339	P1 = [-4.0, -3.0]
1340	P2 = [-3.0, -2.0]
1341	P3 = [-2.0, -1.0]
1342	P4 = [ 3.0, 4.0]
1343	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1344
1345	def test_intersection_bug_20081110_TR_BR(self):
1346	"""test_intersection_bug_20081110(self)
1347
1348	Test all cases in top-right & bottom-right quadrant
1349	"""
1350
1351	# define 4 collinear points, 1 in BR, 3 in TR
1352	# P1-+-P2---P3---P4
1353	P1 = [ 1.0, -3.0]
1354	P2 = [ 5.0, 1.0]
1355	P3 = [ 6.0, 2.0]
1356	P4 = [ 7.0, 3.0]
1357	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1358
1359	# define 4 collinear points, 2 in BR, 2 in TR
1360	# P1---P2-+-P3---P4
1361	P1 = [ 1.0, -3.0]
1362	P2 = [ 2.0, -2.0]
1363	P3 = [ 6.0, 2.0]
1364	P4 = [ 7.0, 3.0]
1365	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1366
1367	# define 4 collinear points, 3 in BR, 1 in TR
1368	# P1---P2---P3-+-P4
1369	P1 = [ 1.0, -3.0]
1370	P2 = [ 2.0, -2.0]
1371	P3 = [ 3.0, -1.0]
1372	P4 = [ 7.0, 3.0]
1373	self.helper_test_parallel_intersection_code(P1, P2, P3, P4)
1374
1375
1376	def test_intersection_direction_invariance(self):
1377	"""This runs through a number of examples and checks that
1378	direction of lines don't matter.
1379	"""
1380
1381	line0 = [[0,0], [100,100]]
1382
1383	common_end_point = [20, 150]
1384
1385	for i in range(100):
1386	x = 20 + i * 1.0/100
1387
1388	line1 = [[x,0], common_end_point]
1389	status, p1 = intersection(line0, line1)
1390	assert status == 1
1391
1392	# Swap direction of line1
1393	line1 = [common_end_point, [x,0]]
1394	status, p2 = intersection(line0, line1)
1395	assert status == 1
1396
1397	msg = ('Orientation of line should not matter.\n'
1398	'However, segment [%f,%f], [%f, %f]' %
1399	(x, 0, common_end_point[0], common_end_point[1]))
1400	msg += ' gave %s, \nbut when reversed we got %s' % (p1, p2)
1401	assert num.allclose(p1, p2), msg
1402
1403	# Swap order of lines
1404	status, p3 = intersection(line1, line0)
1405	assert status == 1
1406	msg = 'Order of lines gave different results'
1407	assert num.allclose(p1, p3), msg
1408
1409	def test_no_intersection(self):
1410	""" Test 2 non-touching lines don't intersect. """
1411	line0 = [[-1,1], [1,1]]
1412	line1 = [[0,-1], [0,0]]
1413
1414	status, value = intersection(line0, line1)
1415	assert status == 0
1416	assert value is None
1417
1418	def test_intersection_parallel(self):
1419	line0 = [[-1,1], [1,1]]
1420	line1 = [[-1,0], [5,0]]
1421
1422	status, value = intersection(line0, line1)
1423	assert status == 4
1424	assert value is None
1425
1426	line0 = [[0,0], [10,100]]
1427	line1 = [[-10,5], [0,105]]
1428
1429	status, value = intersection(line0, line1)
1430	assert status == 4
1431	assert value is None
1432
1433	def test_intersection_coincide(self):
1434	"""Test what happens when two lines partly coincide"""
1435
1436	# Overlap 1
1437	line0 = [[0,0], [5,0]]
1438	line1 = [[-3,0], [3,0]]
1439
1440	status, value = intersection(line0, line1)
1441	assert status == 2
1442	assert num.allclose(value, [[0,0], [3,0]])
1443
1444	# Overlap 2
1445	line0 = [[-10,0], [5,0]]
1446	line1 = [[-3,0], [10,0]]
1447
1448	status, value = intersection(line0, line1)
1449	assert status == 2
1450	assert num.allclose(value, [[-3, 0], [5,0]])
1451
1452	# Inclusion 1
1453	line0 = [[0,0], [5,0]]
1454	line1 = [[2,0], [3,0]]
1455
1456	status, value = intersection(line0, line1)
1457	assert status == 2
1458	assert num.allclose(value, line1)
1459
1460	# Inclusion 2
1461	line0 = [[1,0], [5,0]]
1462	line1 = [[-10,0], [15,0]]
1463
1464	status, value = intersection(line0, line1)
1465	assert status == 2
1466	assert num.allclose(value, line0)
1467
1468	# Exclusion (no intersection)
1469	line0 = [[-10,0], [1,0]]
1470	line1 = [[3,0], [15,0]]
1471
1472	status, value = intersection(line0, line1)
1473	assert status == 3
1474	assert value is None
1475
1476	# Try examples with some slope (y=2*x+5)
1477
1478	# Overlap
1479	line0 = [[0, 5], [7, 19]]
1480	line1 = [[1, 7], [10, 25]]
1481	status, value = intersection(line0, line1)
1482	assert status == 2
1483	assert num.allclose(value, [[1, 7], [7, 19]])
1484
1485	status, value = intersection(line1, line0)
1486	assert status == 2
1487	assert num.allclose(value, [[1, 7], [7, 19]])
1488
1489	# Swap direction
1490	line0 = [[7,19], [0,5]]
1491	line1 = [[1,7], [10,25]]
1492	status, value = intersection(line0, line1)
1493	assert status == 2
1494	assert num.allclose(value, [[7, 19], [1, 7]])
1495
1496	line0 = [[0,5], [7,19]]
1497	line1 = [[10,25], [1,7]]
1498	status, value = intersection(line0, line1)
1499	assert status == 2
1500	assert num.allclose(value, [[1, 7], [7, 19]])
1501
1502	# Inclusion
1503	line0 = [[1,7], [7,19]]
1504	line1 = [[0,5], [10,25]]
1505	status, value = intersection(line0, line1)
1506	assert status == 2
1507	assert num.allclose(value, [[1,7], [7, 19]])
1508
1509	line0 = [[0,5], [10,25]]
1510	line1 = [[1,7], [7,19]]
1511	status, value = intersection(line0, line1)
1512	assert status == 2
1513	assert num.allclose(value, [[1,7], [7, 19]])
1514
1515	line0 = [[0,5], [10,25]]
1516	line1 = [[7,19], [1,7]]
1517	status, value = intersection(line0, line1)
1518	assert status == 2
1519	assert num.allclose(value, [[7, 19], [1, 7]])
1520
1521
1522	def test_inside_polygon_main(self):
1523	# From real example (that failed)
1524	polygon = [[20,20], [40,20], [40,40], [20,40]]
1525	points = [[40, 50]]
1526	res = inside_polygon(points, polygon)
1527	assert len(res) == 0
1528
1529	polygon = [[20,20], [40,20], [40,40], [20,40]]
1530	points = [[25, 25], [30, 20], [40, 50], [90, 20], [40, 90]]
1531	res = inside_polygon(points, polygon)
1532	assert len(res) == 2
1533	assert num.allclose(res, [0,1])
1534
1535	def test_polygon_area(self):
1536	""" Test getting the area of a polygon. """
1537	# Simplest case: Polygon is the unit square
1538	polygon = [[0,0], [1,0], [1,1], [0,1]]
1539	assert polygon_area(polygon) == 1
1540
1541	# Simple case: Polygon is a rectangle
1542	polygon = [[0,0], [1,0], [1,4], [0,4]]
1543	assert polygon_area(polygon) == 4
1544
1545	# Simple case: Polygon is a unit triangle
1546	polygon = [[0,0], [1,0], [0,1]]
1547	assert polygon_area(polygon) == 0.5
1548
1549	# Simple case: Polygon is a diamond
1550	polygon = [[0,0], [1,1], [2,0], [1, -1]]
1551	assert polygon_area(polygon) == 2.0
1552
1553	# Complex case where numerical errors might occur
1554	polygon = [[314037.58727982, 6224952.2960092],
1555	[314038.58727982, 6224952.2960092],
1556	[314038.58727982, 6224953.2960092],
1557	[314037.58727982, 6224953.2960092]]
1558	assert polygon_area(polygon) == 1.0
1559
1560	def test_poly_xy(self):
1561	# Simplest case: Polygon is the unit square
1562	polygon = [[0,0], [1,0], [1,1], [0,1]]
1563	x, y = _poly_xy(polygon)
1564	assert len(x) == len(polygon)+1
1565	assert len(y) == len(polygon)+1
1566	assert x[0] == 0
1567	assert x[1] == 1
1568	assert x[2] == 1
1569	assert x[3] == 0
1570	assert y[0] == 0
1571	assert y[1] == 0
1572	assert y[2] == 1
1573	assert y[3] == 1
1574
1575	# Arbitrary polygon
1576	polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]]
1577	x, y = _poly_xy(polygon)
1578	assert len(x) == len(polygon)+1
1579	assert len(y) == len(polygon)+1
1580	assert x[0] == 1
1581	assert x[1] == 1
1582	assert x[2] == 100
1583	assert x[3] == 1
1584	assert x[4] == 3
1585	assert y[0] == 5
1586	assert y[1] == 1
1587	assert y[2] == 10
1588	assert y[3] == 10
1589	assert y[4] == 6
1590
1591
1592	def test_plot_polygons(self):
1593	import os
1594
1595	# Simplest case: Polygon is the unit square
1596	polygon1 = [[0,0], [1,0], [1,1], [0,1]]
1597	polygon2 = [[1,1], [2,1], [3,2], [2,2]]
1598	plot_polygons([polygon1, polygon2], figname='test1')
1599
1600	# Another case
1601	polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
1602	plot_polygons([polygon2, polygon3], figname='test2')
1603
1604	for file in ['test1.png', 'test2.png']:
1605	assert os.access(file, os.R_OK)
1606	os.remove(file)
1607
1608
1609	def test_inside_polygon_geospatial(self):
1610	""" Test geospatial coords inside poly. """
1611	#Simplest case: Polygon is the unit square
1612	polygon_absolute = [[0, 0], [1, 0], [1, 1], [0, 1]]
1613	poly_geo_ref = Geo_reference(57, 100, 100)
1614	polygon = poly_geo_ref.change_points_geo_ref(polygon_absolute)
1615	poly_spatial = Geospatial_data(polygon, geo_reference=poly_geo_ref)
1616
1617	points_absolute = (0.5, 0.5)
1618	points_geo_ref = Geo_reference(57, 78, -56)
1619	points = points_geo_ref.change_points_geo_ref(points_absolute)
1620	points_spatial = Geospatial_data(points, geo_reference=points_geo_ref)
1621
1622	assert is_inside_polygon(points_absolute, polygon_absolute)
1623	assert is_inside_polygon(ensure_numeric(points_absolute),
1624	ensure_numeric(polygon_absolute))
1625	assert is_inside_polygon(points_absolute, poly_spatial)
1626	assert is_inside_polygon(points_spatial, poly_spatial)
1627	assert is_inside_polygon(points_spatial, polygon_absolute)
1628
1629	assert is_inside_polygon(points_absolute, polygon_absolute)
1630
1631	def test_decimate_polygon(self):
1632	from polygon import decimate_polygon
1633	polygon = [[0,0], [10,10], [15,5], [20, 10],
1634	[25,0], [30,10], [40,-10], [35, -5]]
1635
1636	dpoly = decimate_polygon(polygon, factor=2)
1637
1638	assert len(dpoly)*2==len(polygon)
1639
1640
1641	def test_interpolate_polyline(self):
1642	"""test_interpolate_polyline(self):
1643
1644	This test is added under the assumption that the function
1645	interpolate_polyline implemented by John Jakeman works.
1646	It has been exercised somewhat by tests of sts boundary,
1647	but never before separately.
1648	"""
1649
1650	f = num.array([58.06150614, 58.06150614, 58.06150614])
1651	vertex_coordinates = num.array([[0., 0., ],
1652	[4.04092634, 1106.11074699],
1653	[8.08836552, 2212.16910609]])
1654	gauge_neighbour_id = [1, 2, -1]
1655	point_coordinates = num.array([[2.21870766e+03, 1.09802864e+03],
1656	[1.62739645e+03, 2.20626983e+03],
1657	[5.20084967e+02, 2.21030386e+03],
1658	[6.06464546e+00, 1.65913993e+03],
1659	[1.61934862e+03, -5.88143836e+00],
1660	[5.11996623e+02, -1.85956061e+00],
1661	[2.02046270e+00, 5.53055373e+02]])
1662
1663	z_ref = [0., 0., 0., 58.06150614, 0., 0., 58.06150614]
1664
1665	z = interpolate_polyline(f, vertex_coordinates,
1666	gauge_neighbour_id, point_coordinates)
1667	assert num.allclose(z, z_ref)
1668
1669	# Another f
1670	f = num.array([58.06150614, 158.06150614, 258.06150614])
1671	z_ref = [0., 0., 0., 208.06150645, 0., 0., 108.0615061]
1672	z = interpolate_polyline(f, vertex_coordinates,
1673	gauge_neighbour_id, point_coordinates)
1674	assert num.allclose(z, z_ref)
1675
1676	# Other and simpler numbers
1677	f = num.array([1, 5, 13])
1678	vertex_coordinates = num.array([[0., 0., ],
1679	[4., 4.],
1680	[8., 8.]])
1681	point_coordinates = num.array([[0.1, 0.1],
1682	[3.5, 3.5],
1683	[4.0, 4.0],
1684	[5.2, 5.2],
1685	[7.0, 7.0],
1686	[8.3, 8.3]])
1687	gauge_neighbour_id = [1, 2, -1]
1688
1689	z = interpolate_polyline(f, vertex_coordinates,
1690	gauge_neighbour_id, point_coordinates)
1691	z_ref = [1.1, 4.5, 5., 7.4, 11., 0.]
1692	assert num.allclose(z, z_ref)
1693
1694	# Test exception thrown for one point
1695	f = num.array([5])
1696	vertex_coordinates = num.array([[4., 4.]])
1697	try:
1698	z = interpolate_polyline(f, vertex_coordinates,
1699	gauge_neighbour_id, point_coordinates)
1700	except Exception:
1701	pass
1702	else:
1703	raise Exception, 'One point should have raised exception'
1704
1705	# More polyline nodes
1706	data = num.array([1, 5, 13, 12, 6, 29])
1707	polyline_nodes = num.array([[0., 0.],
1708	[4., 4.],
1709	[8., 8.],
1710	[10., 10.],
1711	[10., 5.],
1712	[10., 0.]])
1713	point_coordinates = num.array([[0.1, 0.1],
1714	[3.5, 3.5],
1715	[4.0, 4.0],
1716	[5.2, 5.2],
1717	[7.0, 7.0],
1718	[8.3, 8.3],
1719	[10., 10.],
1720	[10., 9.],
1721	[10., 7.1],
1722	[10., 4.3],
1723	[10., 1.0]])
1724	gauge_neighbour_id = [1, 2, 3, 4, 5, -1]
1725	z = interpolate_polyline(data, polyline_nodes,
1726	gauge_neighbour_id, point_coordinates)
1727	z_ref = [1.1, 4.5, 5., 7.4, 11., 12.85, 12., 10.8, 8.52, 9.22, 24.4]
1728	assert num.allclose(z, z_ref)
1729
1730
1731	def test_is_inside_triangle_more(self):
1732	""" Test if points inside triangles are detected correctly. """
1733	res = is_inside_triangle([0.5, 0.5], [[ 0.5, 0. ],
1734	[ 0.5, 0.5],
1735	[ 0., 0. ]])
1736	assert res is True
1737
1738	res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1739	[[ 0.5, 0. ], [ 0.5, 0.5], [0., 0.]])
1740	assert res is False
1741
1742	res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1743	[[1., 0.], [1., 0.5], [0.5, 0.]])
1744	assert res is False
1745
1746
1747	res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1748	[[0.5, 0.5], [0.5, 0.], [1., 0.5]])
1749	assert res is True
1750
1751
1752	res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
1753	[[0.5, 0.], [0.5, 0.5], [0., 0.]])
1754	assert res is False
1755
1756
1757	res = is_inside_triangle([0.10000000000000001, 0.20000000000000001],
1758	[[0., 0.5], [0., 0.], [0.5, 0.5]])
1759	assert res is True
1760
1761	res = is_inside_triangle([0.69999999999999996, 0.69999999999999996],
1762	[[0.5, 0.], [0.5, 0.5], [0., 0.]])
1763	assert res is False
1764
1765	res = is_inside_triangle([0.59999999999999998, 0.29999999999999999],
1766	[[0.25, 0.5], [0.25, 0.25], [0.5, 0.5]])
1767	assert res is False
1768
1769
1770	res = is_inside_triangle([10, 3],
1771	[[0.1, 0.],
1772	[0.1, 0.08333333],
1773	[0., 0.]])
1774	assert res is False
1775
1776
1777	def test_is_polygon_complex(self):
1778	""" Test a concave and a complex poly with is_complex, to make
1779	sure it can detect self-intersection.
1780	"""
1781	concave_poly = [[0, 0], [10, 0], [5, 5], [10, 10], [0, 10]]
1782	complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], \
1783	[0, 10]]
1784
1785	assert not is_complex(concave_poly)
1786	assert is_complex(complex_poly)
1787
1788	def test_is_polygon_complex2(self):
1789	""" Test a concave and a complex poly with is_complex, to make sure it
1790	can detect self-intersection. This test uses more complicated
1791	polygons.
1792	"""
1793	concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7, 6], [10, 10], \
1794	[1, 5], [0, 10]]
1795	complex_poly = [[0, 0], [12, 12], [10, 0], [5, 5], [3,18], [4, 15], \
1796	[5, 7], [10, 10], [0, 10], [16, 12]]
1797
1798	assert not is_complex(concave_poly)
1799	assert is_complex(complex_poly)
1800
1801	################################################################################
1802
1803	if __name__ == "__main__":
1804	suite = unittest.makeSuite(Test_Polygon,'test')
1805	runner = unittest.TextTestRunner()
1806	runner.run(suite)

Note: See TracBrowser for help on using the repository browser.

Download in other formats: