Context Navigation

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

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

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

Download in other formats: