Context Navigation

source: branches/numpy/anuga/utilities/test_polygon.py @ 6360

Last change on this file since 6360 was 6360, checked in by rwilson, 15 years ago
Ongoing conversion changes.
File size: 63.0 KB

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

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

Download in other formats: