Context Navigation

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

Last change on this file since 7690 was 7690, checked in by hudson, 15 years ago
Check for complex polygons raises an exception.
File size: 67.7 KB

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

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

Download in other formats: