Context Navigation

test_least_squares.py @ 3236

Last change on this file since 3236 was 3085, checked in by ole, 19 years ago
Moved pyvolution.mesh.py to pyvolution.test_mesh.py in order to avoid confusion with pmesh.
File size: 58.8 KB

Rev	Line
[2526]	1	#!/usr/bin/env python
	2
	3	#TEST
	4	import sys
	5	import unittest
	6	from math import sqrt
	7
	8
[2753]	9	from pyvolution.least_squares import *
[3085]	10	from pyvolution.neighbour_mesh import Mesh
	11
[2526]	12	from Numeric import allclose, array, transpose
	13
	14	from coordinate_transforms.geo_reference import Geo_reference
	15
	16	def distance(x, y):
	17	return sqrt( sum( (array(x)-array(y))**2 ))
	18
	19	def linear_function(point):
	20	point = array(point)
	21	return point[:,0]+point[:,1]
	22
	23
	24	class Test_Least_Squares(unittest.TestCase):
	25
	26	def setUp(self):
	27	pass
	28
	29	def tearDown(self):
	30	pass
	31
	32	def test_datapoint_at_centroid(self):
	33	a = [0.0, 0.0]
	34	b = [0.0, 2.0]
	35	c = [2.0,0.0]
	36	points = [a, b, c]
	37	vertices = [ [1,0,2] ] #bac
	38
	39	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
	40
[2778]	41
[2526]	42	interp = Interpolation(points, vertices, data)
	43	assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]])
	44
	45
	46	def test_quad_tree(self):
	47	p0 = [-10.0, -10.0]
	48	p1 = [20.0, -10.0]
	49	p2 = [-10.0, 20.0]
	50	p3 = [10.0, 50.0]
	51	p4 = [30.0, 30.0]
	52	p5 = [50.0, 10.0]
	53	p6 = [40.0, 60.0]
	54	p7 = [60.0, 40.0]
	55	p8 = [-66.0, 20.0]
	56	p9 = [10.0, -66.0]
	57
	58	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	59	triangles = [ [0, 1, 2],
	60	[3, 2, 4],
	61	[4, 2, 1],
	62	[4, 1, 5],
	63	[3, 4, 6],
	64	[6, 4, 7],
	65	[7, 4, 5],
	66	[8, 0, 2],
	67	[0, 9, 1]]
	68
	69	data = [ [4,4] ]
	70	interp = Interpolation(points, triangles, data, alpha = 0.0,
	71	max_points_per_cell = 4)
	72	#print "PDSG - interp.get_A()", interp.get_A()
	73	answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0.,
	74	0., 0. , 0., 0., 0., 0.]]
	75	assert allclose(interp.get_A(), answer)
	76	interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	77	#print "PDSG - interp.get_A()", interp.get_A()
	78	answer = [ [ 0.0, 0.0, 0.0, 0.,
	79	0., 0. , 0., 0., 0., 0.]]
	80	assert allclose(interp.get_A(), answer)
	81
	82
	83	#point outside of quad tree root cell
	84	interp.set_point_coordinates([[-70, -70]])
	85	#print "PDSG - interp.get_A()", interp.get_A()
	86	answer = [ [ 0.0, 0.0, 0.0, 0.,
	87	0., 0. , 0., 0., 0., 0.]]
	88	assert allclose(interp.get_A(), answer)
	89
	90	def test_expand_search(self):
	91	p0 = [-10.0, -10.0]
	92	p1 = [20.0, -10.0]
	93	p2 = [-10.0, 20.0]
	94	p3 = [10.0, 50.0]
	95	p4 = [30.0, 30.0]
	96	p5 = [50.0, 10.0]
	97	p6 = [40.0, 60.0]
	98	p7 = [60.0, 40.0]
	99	p8 = [-66.0, 20.0]
	100	p9 = [10.0, -66.0]
	101
	102	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	103	triangles = [ [0, 1, 2],
	104	[3, 2, 4],
	105	[4, 2, 1],
	106	[4, 1, 5],
	107	[3, 4, 6],
	108	[6, 4, 7],
	109	[7, 4, 5],
	110	[8, 0, 2],
	111	[0, 9, 1]]
	112
	113	data = [ [4,4],
	114	[-30,10],
	115	[-20,0],
	116	[-20,10],
	117	[0,30],
	118	[10,-40],
	119	[10,-30],
	120	[10,-20],
	121	[10,10],
	122	[10,20],
	123	[10,30],
	124	[10,40],
	125	[20,10],
	126	[25,45],
	127	[30,0],
	128	[30,10],
	129	[30,30],
	130	[30,40],
	131	[30,50],
	132	[40,10],
	133	[40,30],
	134	[40,40],
	135	[40,50],
	136	[50,20],
	137	[50,30],
	138	[50,40],
	139	[50,50],
	140	[30,0],
	141	[-20,-20]]
	142	point_attributes = [ -400000,
	143	10,
	144	10,
	145	10,
	146	10,
	147	10,
	148	10,
	149	10,
	150	10,
	151	10,
	152	10,
	153	10,
	154	10,
	155	10,
	156	10,
	157	10,
	158	10,
	159	10,
	160	10,
	161	10,
	162	10,
	163	10,
	164	10,
	165	10,
	166	10,
	167	10,
	168	10,
	169	10,
	170	99]
	171
	172	interp = Interpolation(points, triangles, data,
	173	alpha=0.0, expand_search=False, #verbose = True, #False,
	174	max_points_per_cell = 4)
	175	calc = interp.fit_points(point_attributes, )
	176	#print "calc",calc
	177
	178	# the point at 4,4 is ignored. An expanded search has to be done
	179	# to fine which triangel it's in.
	180	# An expanded search isn't done to find that the last point
	181	# isn't in the mesh. But this isn't tested.
	182	answer= [ 10,
	183	10,
	184	10,
	185	10,
	186	10,
	187	10,
	188	10,
	189	10,
	190	10,
	191	10]
	192	assert allclose(calc, answer)
	193
	194	def test_quad_treeII(self):
	195	p0 = [-66.0, 14.0]
	196	p1 = [14.0, -66.0]
	197	p2 = [14.0, 14.0]
	198	p3 = [60.0, 20.0]
	199	p4 = [10.0, 60.0]
	200	p5 = [60.0, 60.0]
	201
	202	points = [p0, p1, p2, p3, p4, p5]
	203	triangles = [ [0, 1, 2],
	204	[3, 2, 1],
	205	[0, 2, 4],
	206	[4, 2, 5],
	207	[5, 2, 3]]
	208
	209	data = [ [-26.0,-26.0] ]
	210	interp = Interpolation(points, triangles, data, alpha = 0.0,
	211	max_points_per_cell = 4)
	212	#print "PDSG - interp.get_A()", interp.get_A()
	213	answer = [ [ 0.5, 0.5, 0.0, 0.,
	214	0., 0.]]
	215	assert allclose(interp.get_A(), answer)
	216	interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	217	#print "PDSG -30,-30 - interp.get_A()", interp.get_A()
	218	answer = [ [ 0.0, 0.0, 0.0, 0.,
	219	0., 0.]]
	220	assert allclose(interp.get_A(), answer)
	221
	222
	223	#point outside of quad tree root cell
	224	interp.set_point_coordinates([[-70, -70]])
	225	#print "PDSG -70,-70 interp.get_A()", interp.get_A()
	226	answer = [ [ 0.0, 0.0, 0.0, 0.,
	227	0., 0. ]]
	228	assert allclose(interp.get_A(), answer)
	229
	230
	231	def test_datapoints_at_vertices(self):
	232	"""Test that data points coinciding with vertices yield a diagonal matrix
	233	"""
	234
	235	a = [0.0, 0.0]
	236	b = [0.0, 2.0]
	237	c = [2.0,0.0]
	238	points = [a, b, c]
	239	vertices = [ [1,0,2] ] #bac
	240
	241	data = points #Use data at vertices
	242
	243	interp = Interpolation(points, vertices, data)
	244	assert allclose(interp.get_A(), [[1., 0., 0.],
	245	[0., 1., 0.],
	246	[0., 0., 1.]])
	247
	248
	249
	250	def test_datapoints_on_edge_midpoints(self):
	251	"""Try datapoints midway on edges -
	252	each point should affect two matrix entries equally
	253	"""
	254
	255	a = [0.0, 0.0]
	256	b = [0.0, 2.0]
	257	c = [2.0,0.0]
	258	points = [a, b, c]
	259	vertices = [ [1,0,2] ] #bac
	260
	261	data = [ [0., 1.], [1., 0.], [1., 1.] ]
	262
	263	interp = Interpolation(points, vertices, data)
	264
	265	assert allclose(interp.get_A(), [[0.5, 0.5, 0.0], #Affects vertex 1 and 0
	266	[0.5, 0.0, 0.5], #Affects vertex 0 and 2
	267	[0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
	268
	269
	270	def test_datapoints_on_edges(self):
	271	"""Try datapoints on edges -
	272	each point should affect two matrix entries in proportion
	273	"""
	274
	275	a = [0.0, 0.0]
	276	b = [0.0, 2.0]
	277	c = [2.0,0.0]
	278	points = [a, b, c]
	279	vertices = [ [1,0,2] ] #bac
	280
	281	data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
	282
	283	interp = Interpolation(points, vertices, data)
	284
	285	assert allclose(interp.get_A(), [[0.25, 0.75, 0.0], #Affects vertex 1 and 0
	286	[0.25, 0.0, 0.75], #Affects vertex 0 and 2
	287	[0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
	288
	289	def test_arbitrary_datapoints(self):
	290	"""Try arbitrary datapoints
	291	"""
	292
	293	from Numeric import sum
	294
	295	a = [0.0, 0.0]
	296	b = [0.0, 2.0]
	297	c = [2.0,0.0]
	298	points = [a, b, c]
	299	vertices = [ [1,0,2] ] #bac
	300
	301	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
	302
	303	interp = Interpolation(points, vertices, data)
	304	#print "interp.get_A()", interp.get_A()
	305	assert allclose(sum(interp.get_A(), axis=1), 1.0)
	306
	307	def test_arbitrary_datapoints_some_outside(self):
	308	"""Try arbitrary datapoints one outside the triangle.
	309	That one should be ignored
	310	"""
	311
	312	from Numeric import sum
	313
	314	a = [0.0, 0.0]
	315	b = [0.0, 2.0]
	316	c = [2.0,0.0]
	317	points = [a, b, c]
	318	vertices = [ [1,0,2] ] #bac
	319
	320	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
	321
	322
	323	interp = Interpolation(points, vertices, data, precrop = True)
	324	assert allclose(sum(interp.get_A(), axis=1), 1.0)
	325
	326	interp = Interpolation(points, vertices, data, precrop = False)
	327	assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0])
	328
	329
	330
	331	# this causes a memory error in scipy.sparse
	332	def test_more_triangles(self):
	333
	334	a = [-1.0, 0.0]
	335	b = [3.0, 4.0]
	336	c = [4.0,1.0]
	337	d = [-3.0, 2.0] #3
	338	e = [-1.0,-2.0]
	339	f = [1.0, -2.0] #5
	340
	341	points = [a, b, c, d,e,f]
	342	triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
	343
	344	#Data points
	345	data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
	346	interp = Interpolation(points, triangles, data_points)
	347
	348	answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d
	349	[0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	350	[0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	351	[0.0, 0.5, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	352	[0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	353	[1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
	354
	355
	356	A = interp.get_A()
	357	for i in range(A.shape[0]):
	358	for j in range(A.shape[1]):
	359	if not allclose(A[i,j], answer[i][j]):
	360	print i,j,':',A[i,j], answer[i][j]
	361
	362
	363	assert allclose(interp.get_A(), answer)
	364
	365
	366
	367
	368	def test_smooth_attributes_to_mesh(self):
	369	a = [0.0, 0.0]
	370	b = [0.0, 5.0]
	371	c = [5.0, 0.0]
	372	points = [a, b, c]
	373	triangles = [ [1,0,2] ] #bac
	374
	375	d1 = [1.0, 1.0]
	376	d2 = [1.0, 3.0]
	377	d3 = [3.0,1.0]
	378	z1 = 2
	379	z2 = 4
	380	z3 = 4
	381	data_coords = [d1, d2, d3]
	382
	383	interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
	384	z = [z1, z2, z3]
	385	f = interp.fit(z)
	386	answer = [0, 5., 5.]
	387
	388	#print "f\n",f
	389	#print "answer\n",answer
	390
	391	assert allclose(f, answer, atol=1e-7)
	392
	393
	394	def test_smooth_att_to_meshII(self):
	395
	396	a = [0.0, 0.0]
	397	b = [0.0, 5.0]
	398	c = [5.0, 0.0]
	399	points = [a, b, c]
	400	triangles = [ [1,0,2] ] #bac
	401
	402	d1 = [1.0, 1.0]
	403	d2 = [1.0, 2.0]
	404	d3 = [3.0,1.0]
	405	data_coords = [d1, d2, d3]
	406	z = linear_function(data_coords)
	407	interp = Interpolation(points, triangles, data_coords, alpha=0.0)
	408	f = interp.fit(z)
	409	answer = linear_function(points)
	410
	411	assert allclose(f, answer)
	412
	413	def test_smooth_attributes_to_meshIII(self):
	414
	415	a = [-1.0, 0.0]
	416	b = [3.0, 4.0]
	417	c = [4.0,1.0]
	418	d = [-3.0, 2.0] #3
	419	e = [-1.0,-2.0]
	420	f = [1.0, -2.0] #5
	421
	422	vertices = [a, b, c, d,e,f]
	423	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	424
	425	point_coords = [[-2.0, 2.0],
	426	[-1.0, 1.0],
	427	[0.0,2.0],
	428	[1.0, 1.0],
	429	[2.0, 1.0],
	430	[0.0,0.0],
	431	[1.0, 0.0],
	432	[0.0, -1.0],
	433	[-0.2,-0.5],
	434	[-0.9, -1.5],
	435	[0.5, -1.9],
	436	[3.0,1.0]]
	437
	438	z = linear_function(point_coords)
	439	interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
	440
	441	#print 'z',z
	442	f = interp.fit(z)
	443	answer = linear_function(vertices)
	444	#print "f\n",f
	445	#print "answer\n",answer
	446	assert allclose(f, answer)
	447
	448
	449	def test_smooth_attributes_to_meshIV(self):
	450	""" Testing 2 attributes smoothed to the mesh
	451	"""
	452
	453	a = [0.0, 0.0]
	454	b = [0.0, 5.0]
	455	c = [5.0, 0.0]
	456	points = [a, b, c]
	457	triangles = [ [1,0,2] ] #bac
	458
	459	d1 = [1.0, 1.0]
	460	d2 = [1.0, 3.0]
	461	d3 = [3.0, 1.0]
	462	z1 = [2, 4]
	463	z2 = [4, 8]
	464	z3 = [4, 8]
	465	data_coords = [d1, d2, d3]
	466
	467	interp = Interpolation(points, triangles, data_coords, alpha=0.0)
	468	z = [z1, z2, z3]
	469	f = interp.fit_points(z)
	470	answer = [[0,0], [5., 10.], [5., 10.]]
	471	assert allclose(f, answer)
	472
	473	def test_interpolate_attributes_to_points(self):
	474	v0 = [0.0, 0.0]
	475	v1 = [0.0, 5.0]
	476	v2 = [5.0, 0.0]
	477
	478	vertices = [v0, v1, v2]
	479	triangles = [ [1,0,2] ] #bac
	480
	481	d0 = [1.0, 1.0]
	482	d1 = [1.0, 2.0]
	483	d2 = [3.0, 1.0]
	484	point_coords = [ d0, d1, d2]
	485
	486	interp = Interpolation(vertices, triangles, point_coords)
	487	f = linear_function(vertices)
[2750]	488	#z = interp.interpolate(f)
[2526]	489	answer = linear_function(point_coords)
	490
	491
[2750]	492	#assert allclose(z, answer)
[2526]	493
	494
	495	def test_interpolate_attributes_to_points_interp_only(self):
	496	v0 = [0.0, 0.0]
	497	v1 = [0.0, 5.0]
	498	v2 = [5.0, 0.0]
	499
	500	vertices = [v0, v1, v2]
	501	triangles = [ [1,0,2] ] #bac
	502
	503	d0 = [1.0, 1.0]
	504	d1 = [1.0, 2.0]
	505	d2 = [3.0, 1.0]
	506	point_coords = [ d0, d1, d2]
	507
	508	interp = Interpolation(vertices, triangles, point_coords,
	509	interp_only = True)
	510
	511	f = linear_function(vertices)
[2750]	512	#z = interp.interpolate(f)
[2526]	513	answer = linear_function(point_coords)
	514	#print "answer", answer
	515	#print "z", z
	516
[2750]	517	#assert allclose(z, answer)
[2526]	518
	519	def test_interpolate_attributes_to_pointsII(self):
	520	a = [-1.0, 0.0]
	521	b = [3.0, 4.0]
	522	c = [4.0, 1.0]
	523	d = [-3.0, 2.0] #3
	524	e = [-1.0, -2.0]
	525	f = [1.0, -2.0] #5
	526
	527	vertices = [a, b, c, d,e,f]
	528	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	529
	530
	531	point_coords = [[-2.0, 2.0],
	532	[-1.0, 1.0],
	533	[0.0, 2.0],
	534	[1.0, 1.0],
	535	[2.0, 1.0],
	536	[0.0, 0.0],
	537	[1.0, 0.0],
	538	[0.0, -1.0],
	539	[-0.2, -0.5],
	540	[-0.9, -1.5],
	541	[0.5, -1.9],
	542	[3.0, 1.0]]
	543
	544	interp = Interpolation(vertices, triangles, point_coords)
	545	f = linear_function(vertices)
[2750]	546	#z = interp.interpolate(f)
[2526]	547	answer = linear_function(point_coords)
	548	#print "z",z
	549	#print "answer",answer
[2750]	550	#assert allclose(z, answer)
[2526]	551
	552	def test_interpolate_attributes_to_pointsIII(self):
	553	"""Test linear interpolation of known values at vertices to
	554	new points inside a triangle
	555	"""
	556	a = [0.0, 0.0]
	557	b = [0.0, 5.0]
	558	c = [5.0, 0.0]
	559	d = [5.0, 5.0]
	560
	561	vertices = [a, b, c, d]
	562	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	563
	564	#Points within triangle 1
	565	d0 = [1.0, 1.0]
	566	d1 = [1.0, 2.0]
	567	d2 = [3.0, 1.0]
	568
	569	#Point within triangle 2
	570	d3 = [4.0, 3.0]
	571
	572	#Points on common edge
	573	d4 = [2.5, 2.5]
	574	d5 = [4.0, 1.0]
	575
	576	#Point on common vertex
	577	d6 = [0., 5.]
	578
	579	point_coords = [d0, d1, d2, d3, d4, d5, d6]
	580
	581	interp = Interpolation(vertices, triangles, point_coords)
	582
	583	#Known values at vertices
	584	#Functions are x+y, x+2y, 2x+y, x-y-5
	585	f = [ [0., 0., 0., -5.], # (0,0)
	586	[5., 10., 5., -10.], # (0,5)
	587	[5., 5., 10.0, 0.], # (5,0)
	588	[10., 15., 15., -5.]] # (5,5)
	589
[2750]	590	#z = interp.interpolate(f)
[2526]	591	answer = [ [2., 3., 3., -5.], # (1,1)
	592	[3., 5., 4., -6.], # (1,2)
	593	[4., 5., 7., -3.], # (3,1)
	594	[7., 10., 11., -4.], # (4,3)
	595	[5., 7.5, 7.5, -5.], # (2.5, 2.5)
	596	[5., 6., 9., -2.], # (4,1)
	597	[5., 10., 5., -10.]] # (0,5)
	598
	599	#print "***********"
	600	#print "z",z
	601	#print "answer",answer
	602	#print "***********"
	603
	604	#Should an error message be returned if points are outside
	605	# of the mesh? Not currently.
	606
[2750]	607	#assert allclose(z, answer)
[2526]	608
	609
	610	def test_interpolate_point_outside_of_mesh(self):
	611	"""Test linear interpolation of known values at vertices to
	612	new points inside a triangle
	613	"""
	614	a = [0.0, 0.0]
	615	b = [0.0, 5.0]
	616	c = [5.0, 0.0]
	617	d = [5.0, 5.0]
	618
	619	vertices = [a, b, c, d]
	620	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	621
	622	#Far away point
	623	d7 = [-1., -1.]
	624
	625	point_coords = [ d7]
	626
	627	interp = Interpolation(vertices, triangles, point_coords)
	628
	629	#Known values at vertices
	630	#Functions are x+y, x+2y, 2x+y, x-y-5
	631	f = [ [0., 0., 0., -5.], # (0,0)
	632	[5., 10., 5., -10.], # (0,5)
	633	[5., 5., 10.0, 0.], # (5,0)
	634	[10., 15., 15., -5.]] # (5,5)
	635
	636
[2750]	637	#z = interp.interpolate(f)
	638	#answer = [ [0., 0., 0., 0.]] # (-1,-1)
	639
[2526]	640	#print "***********"
	641	#print "z",z
	642	#print "answer",answer
	643	#print "***********"
	644
	645	#Should an error message be returned if points are outside
	646	# of the mesh? Not currently.
	647
[2750]	648	#assert allclose(z, answer)
[2526]	649
	650	def test_interpolate_attributes_to_pointsIV(self):
	651	a = [-1.0, 0.0]
	652	b = [3.0, 4.0]
	653	c = [4.0, 1.0]
	654	d = [-3.0, 2.0] #3
	655	e = [-1.0, -2.0]
	656	f = [1.0, -2.0] #5
	657
	658	vertices = [a, b, c, d,e,f]
	659	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	660
	661
	662	point_coords = [[-2.0, 2.0],
	663	[-1.0, 1.0],
	664	[0.0, 2.0],
	665	[1.0, 1.0],
	666	[2.0, 1.0],
	667	[0.0, 0.0],
	668	[1.0, 0.0],
	669	[0.0, -1.0],
	670	[-0.2, -0.5],
	671	[-0.9, -1.5],
	672	[0.5, -1.9],
	673	[3.0, 1.0]]
	674
	675	interp = Interpolation(vertices, triangles, point_coords)
	676	f = array([linear_function(vertices),2*linear_function(vertices) ])
	677	f = transpose(f)
	678	#print "f",f
[2750]	679	#z = interp.interpolate(f)
[2526]	680	answer = [linear_function(point_coords),
	681	2*linear_function(point_coords) ]
	682	answer = transpose(answer)
	683	#print "z",z
	684	#print "answer",answer
[2750]	685	#assert allclose(z, answer)
[2526]	686
	687	def test_smooth_attributes_to_mesh_function(self):
	688	""" Testing 2 attributes smoothed to the mesh
	689	"""
	690
	691	a = [0.0, 0.0]
	692	b = [0.0, 5.0]
	693	c = [5.0, 0.0]
	694	points = [a, b, c]
	695	triangles = [ [1,0,2] ] #bac
	696
	697	d1 = [1.0, 1.0]
	698	d2 = [1.0, 3.0]
	699	d3 = [3.0, 1.0]
	700	z1 = [2, 4]
	701	z2 = [4, 8]
	702	z3 = [4, 8]
	703	data_coords = [d1, d2, d3]
	704	z = [z1, z2, z3]
	705
	706	f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
	707	answer = [[0, 0], [5., 10.], [5., 10.]]
	708
	709	assert allclose(f, answer)
	710
	711
	712
	713	def test_pts2rectangular(self):
	714
	715	import time, os
	716	FN = 'xyatest' + str(time.time()) + '.xya'
	717	fid = open(FN, 'w')
	718	fid.write(' %s \n' %('elevation'))
	719	fid.write('%f %f %f\n' %(1,1,2) )
	720	fid.write('%f %f %f\n' %(1,3,4) )
	721	fid.write('%f %f %f\n' %(3,1,4) )
	722	fid.close()
	723
	724	points, triangles, boundary, attributes =\
	725	pts2rectangular(FN, 4, 4)
	726
	727
	728	data_coords = [ [1,1], [1,3], [3,1] ]
	729	z = [2, 4, 4]
	730
[2585]	731	ref = fit_to_mesh(points, triangles, data_coords, z, verbose=False)
[2526]	732
	733	#print attributes
	734	#print ref
	735	assert allclose(attributes, ref)
	736
	737	os.remove(FN)
	738
	739
	740	#Tests of smoothing matrix
	741	def test_smoothing_matrix_one_triangle(self):
	742	from Numeric import dot
	743	a = [0.0, 0.0]
	744	b = [0.0, 2.0]
	745	c = [2.0,0.0]
	746	points = [a, b, c]
	747
	748	vertices = [ [1,0,2] ] #bac
	749
	750	interp = Interpolation(points, vertices)
	751
	752	assert allclose(interp.get_D(), [[1, -0.5, -0.5],
	753	[-0.5, 0.5, 0],
	754	[-0.5, 0, 0.5]])
	755
	756	#Define f(x,y) = x
	757	f = array([0,0,2]) #Value at global vertex 2
	758
	759	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	760	# int 1 dx dy = area = 2
	761	assert dot(dot(f, interp.get_D()), f) == 2
	762
	763	#Define f(x,y) = y
	764	f = array([0,2,0]) #Value at global vertex 1
	765
	766	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	767	# int 1 dx dy = area = 2
	768	assert dot(dot(f, interp.get_D()), f) == 2
	769
	770	#Define f(x,y) = x+y
	771	f = array([0,2,2]) #Values at global vertex 1 and 2
	772
	773	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	774	# int 2 dx dy = 2*area = 4
	775	assert dot(dot(f, interp.get_D()), f) == 4
	776
	777
	778
	779	def test_smoothing_matrix_more_triangles(self):
	780	from Numeric import dot
	781
	782	a = [0.0, 0.0]
	783	b = [0.0, 2.0]
	784	c = [2.0,0.0]
	785	d = [0.0, 4.0]
	786	e = [2.0, 2.0]
	787	f = [4.0,0.0]
	788
	789	points = [a, b, c, d, e, f]
	790	#bac, bce, ecf, dbe, daf, dae
	791	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	792
	793	interp = Interpolation(points, vertices)
	794
	795
	796	#assert allclose(interp.get_D(), [[1, -0.5, -0.5],
	797	# [-0.5, 0.5, 0],
	798	# [-0.5, 0, 0.5]])
	799
	800	#Define f(x,y) = x
	801	f = array([0,0,2,0,2,4]) #f evaluated at points a-f
	802
	803	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	804	# int 1 dx dy = total area = 8
	805	assert dot(dot(f, interp.get_D()), f) == 8
	806
	807	#Define f(x,y) = y
	808	f = array([0,2,0,4,2,0]) #f evaluated at points a-f
	809
	810	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	811	# int 1 dx dy = area = 8
	812	assert dot(dot(f, interp.get_D()), f) == 8
	813
	814	#Define f(x,y) = x+y
	815	f = array([0,2,2,4,4,4]) #f evaluated at points a-f
	816
	817	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	818	# int 2 dx dy = 2*area = 16
	819	assert dot(dot(f, interp.get_D()), f) == 16
	820
	821
	822	def test_fit_and_interpolation(self):
	823
	824	a = [0.0, 0.0]
	825	b = [0.0, 2.0]
	826	c = [2.0, 0.0]
	827	d = [0.0, 4.0]
	828	e = [2.0, 2.0]
	829	f = [4.0, 0.0]
	830
	831	points = [a, b, c, d, e, f]
	832	#bac, bce, ecf, dbe, daf, dae
	833	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	834
	835	#Get (enough) datapoints
	836	data_points = [[ 0.66666667, 0.66666667],
	837	[ 1.33333333, 1.33333333],
	838	[ 2.66666667, 0.66666667],
	839	[ 0.66666667, 2.66666667],
	840	[ 0.0, 1.0],
	841	[ 0.0, 3.0],
	842	[ 1.0, 0.0],
	843	[ 1.0, 1.0],
	844	[ 1.0, 2.0],
	845	[ 1.0, 3.0],
	846	[ 2.0, 1.0],
	847	[ 3.0, 0.0],
	848	[ 3.0, 1.0]]
	849
	850	interp = Interpolation(points, triangles, data_points, alpha=0.0)
	851
	852	z = linear_function(data_points)
	853	answer = linear_function(points)
	854
	855	f = interp.fit(z)
	856
	857	#print "f",f
	858	#print "answer",answer
	859	assert allclose(f, answer)
	860
	861	#Map back
[2750]	862	#z1 = interp.interpolate(f)
[2526]	863	#print "z1\n", z1
	864	#print "z\n",z
[2750]	865	#assert allclose(z, z1)
[2526]	866
	867
	868	def test_smoothing_and_interpolation(self):
	869
	870	a = [0.0, 0.0]
	871	b = [0.0, 2.0]
	872	c = [2.0, 0.0]
	873	d = [0.0, 4.0]
	874	e = [2.0, 2.0]
	875	f = [4.0, 0.0]
	876
	877	points = [a, b, c, d, e, f]
	878	#bac, bce, ecf, dbe, daf, dae
	879	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	880
	881	#Get (too few!) datapoints
	882	data_points = [[ 0.66666667, 0.66666667],
	883	[ 1.33333333, 1.33333333],
	884	[ 2.66666667, 0.66666667],
	885	[ 0.66666667, 2.66666667]]
	886
	887	z = linear_function(data_points)
	888	answer = linear_function(points)
	889
	890	#Make interpolator with too few data points and no smoothing
	891	interp = Interpolation(points, triangles, data_points, alpha=0.0)
	892	#Must raise an exception
	893	try:
	894	f = interp.fit(z)
	895	except:
	896	pass
	897
	898	#Now try with smoothing parameter
	899	interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
	900
	901	f = interp.fit(z)
	902	#f will be different from answer due to smoothing
	903	assert allclose(f, answer,atol=5)
	904
	905	#Map back
[2750]	906	#z1 = interp.interpolate(f)
	907	#assert allclose(z, z1)
[2526]	908
	909
	910
	911	def test_fit_and_interpolation_with_new_points(self):
	912	"""Fit a surface to one set of points. Then interpolate that surface
	913	using another set of points.
	914	"""
	915
	916	#Setup mesh used to represent fitted function
	917	a = [0.0, 0.0]
	918	b = [0.0, 2.0]
	919	c = [2.0, 0.0]
	920	d = [0.0, 4.0]
	921	e = [2.0, 2.0]
	922	f = [4.0, 0.0]
	923
	924	points = [a, b, c, d, e, f]
	925	#bac, bce, ecf, dbe, daf, dae
	926	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	927
	928	#Datapoints to fit from
	929	data_points1 = [[ 0.66666667, 0.66666667],
	930	[ 1.33333333, 1.33333333],
	931	[ 2.66666667, 0.66666667],
	932	[ 0.66666667, 2.66666667],
	933	[ 0.0, 1.0],
	934	[ 0.0, 3.0],
	935	[ 1.0, 0.0],
	936	[ 1.0, 1.0],
	937	[ 15, -17], #Outside mesh
	938	[ 1.0, 2.0],
	939	[ 1.0, 3.0],
	940	[ 2.0, 1.0],
	941	[ 3.0, 0.0],
	942	[ 3.0, 1.0]]
	943
	944	#Fit surface to mesh
	945	interp = Interpolation(points, triangles, data_points1, alpha=0.0,
[2585]	946	precrop = True, verbose=False)
[2526]	947	z = linear_function(data_points1) #Example z-values
	948	f = interp.fit(z) #Fitted values at vertices
	949
	950
	951
	952	#New datapoints where interpolated values are sought
	953	data_points2 = [[ 0.0, 0.0],
	954	[ 0.5, 0.5],
	955	[ 0.7, 0.7],
	956	[-13, 65], #Outside
	957	[ 1.0, 0.5],
	958	[ 2.0, 0.4],
	959	[ 2.8, 1.2]]
	960
	961
	962
	963	#Build new A matrix based on new points (without precrop)
	964	interp.build_interpolation_matrix_A(data_points2, precrop = False)
	965
	966	#Interpolate using fitted surface
[2750]	967	#z1 = interp.interpolate(f)
[2526]	968
	969	#import Numeric
	970	#data_points2 = Numeric.take(data_points2, interp.point_indices)
	971
	972	#Desired result (OK for points inside)
	973
	974	answer = linear_function(data_points2)
	975	import Numeric
[2750]	976	#z1 = Numeric.take(z1, [0,1,2,4,5,6])
[2526]	977	answer = Numeric.take(answer, [0,1,2,4,5,6])
[2750]	978	#assert allclose(z1, answer)
[2526]	979
	980	#Build new A matrix based on new points (with precrop)
	981	interp.build_interpolation_matrix_A(data_points2, precrop = True)
	982
	983	#Interpolate using fitted surface
[2750]	984	#z1 = interp.interpolate(f)
[2526]	985
	986	import Numeric
	987	data_points2 = Numeric.take(data_points2, interp.point_indices)
	988
	989	#Desired result
	990	answer = linear_function(data_points2)
[2750]	991	#assert allclose(z1, answer)
[2526]	992
	993
	994
	995
	996
	997
	998	def test_interpolation_from_discontinuous_vertex_values(self):
	999	"""test_interpolation_from_discontinuous_vertex_values.
	1000	This will test the format used internally in pyvolution and also
	1001	interpolation from sww files
	1002	"""
	1003
	1004	#Setup mesh used to represent discontinuous function
	1005	a = [0.0, 0.0]
	1006	b = [0.0, 2.0]
	1007	c = [2.0, 0.0]
	1008	d = [0.0, 4.0]
	1009	e = [2.0, 2.0]
	1010	f = [4.0, 0.0]
	1011
	1012	points = [b, a, c,
	1013	b, c, e,
	1014	e, c, f,
	1015	d, b, e]
	1016
	1017	#bac, bce, ecf, dbe
	1018	triangles = [[0,1,2], [3,4,5], [6,7,8], [9,10,11]]
	1019
	1020
	1021	vertex_values = [0.,0.,0.,1.,1.,1.,2.,2.,2.,7.,3.,3.]
	1022
	1023
	1024
	1025	#New datapoints where interpolated values are sought
	1026	data_points = [[0.0, 0.0], #T0
	1027	[0.5, 0.5], #T0
	1028	[1.5, 1.5], #T1
	1029	[2.5, 0.5], #T2
	1030	[0.0, 3.0], #T3
	1031	[1.0, 2.0], #In between T1 and T3 (T1 is used) FIXME?
	1032	[2.0, 1.0], #In between T1 and T2 (T1 is used) FIXME?
	1033	[1.0, 1.0]] #In between T1 and T0 (T0 is used) FIXME?
	1034
	1035
	1036
	1037
	1038	#Build interpolation matrix
	1039	interp = Interpolation(points, triangles, data_points)
	1040	#, alpha=0.0, precrop = True)
	1041
	1042	#print interp.A.todense()
	1043	#print vertex_values
	1044
	1045	#Interpolate using fitted surface
[2750]	1046	#z = interp.interpolate(vertex_values)
[2526]	1047
	1048	#print z
	1049
[2750]	1050	#assert allclose(z, [0,0,1,2,5,1,1,0])
[2526]	1051
	1052
	1053
	1054
	1055	def test_interpolation_function_time_only(self):
	1056	"""Test spatio-temporal interpolation
	1057	Test that spatio temporal function performs the correct
	1058	interpolations in both time and space
	1059	"""
	1060
	1061
	1062	#Three timesteps
	1063	time = [1.0, 5.0, 6.0]
	1064
	1065
	1066	#One quantity
	1067	Q = zeros( (3,6), Float )
	1068
	1069	#Linear in time and space
	1070	a = [0.0, 0.0]
	1071	b = [0.0, 2.0]
	1072	c = [2.0, 0.0]
	1073	d = [0.0, 4.0]
	1074	e = [2.0, 2.0]
	1075	f = [4.0, 0.0]
	1076
	1077	points = [a, b, c, d, e, f]
	1078
	1079	for i, t in enumerate(time):
	1080	Q[i, :] = t*linear_function(points)
	1081
	1082
	1083	#Check basic interpolation of one quantity using averaging
	1084	#(no interpolation points or spatial info)
	1085	from utilities.numerical_tools import mean
	1086	I = Interpolation_function(time, [mean(Q[0,:]),
	1087	mean(Q[1,:]),
	1088	mean(Q[2,:])])
	1089
	1090
	1091
	1092	#Check temporal interpolation
	1093	for i in [0,1,2]:
	1094	assert allclose(I(time[i]), mean(Q[i,:]))
	1095
	1096	#Midway
	1097	assert allclose(I( (time[0] + time[1])/2 ),
	1098	(I(time[0]) + I(time[1]))/2 )
	1099
	1100	assert allclose(I( (time[1] + time[2])/2 ),
	1101	(I(time[1]) + I(time[2]))/2 )
	1102
	1103	assert allclose(I( (time[0] + time[2])/2 ),
	1104	(I(time[0]) + I(time[2]))/2 )
	1105
	1106	#1/3
	1107	assert allclose(I( (time[0] + time[2])/3 ),
	1108	(I(time[0]) + I(time[2]))/3 )
	1109
	1110
	1111	#Out of bounds checks
	1112	try:
	1113	I(time[0]-1)
	1114	except:
	1115	pass
	1116	else:
	1117	raise 'Should raise exception'
	1118
	1119	try:
	1120	I(time[-1]+1)
	1121	except:
	1122	pass
	1123	else:
	1124	raise 'Should raise exception'
	1125
	1126
	1127
	1128
[2750]	1129	def interpolation_test_interpolation_function_spatial_only(self):
[2526]	1130	"""Test spatio-temporal interpolation with constant time
	1131	"""
	1132
	1133	#Three timesteps
	1134	time = [1.0, 5.0, 6.0]
	1135
	1136
	1137	#Setup mesh used to represent fitted function
	1138	a = [0.0, 0.0]
	1139	b = [0.0, 2.0]
	1140	c = [2.0, 0.0]
	1141	d = [0.0, 4.0]
	1142	e = [2.0, 2.0]
	1143	f = [4.0, 0.0]
	1144
	1145	points = [a, b, c, d, e, f]
	1146	#bac, bce, ecf, dbe
	1147	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1148
	1149
	1150	#New datapoints where interpolated values are sought
	1151	interpolation_points = [[ 0.0, 0.0],
	1152	[ 0.5, 0.5],
	1153	[ 0.7, 0.7],
	1154	[ 1.0, 0.5],
	1155	[ 2.0, 0.4],
	1156	[ 2.8, 1.2]]
	1157
	1158
	1159	#One quantity linear in space
	1160	Q = linear_function(points)
	1161
	1162
	1163	#Check interpolation of one quantity using interpolaton points
	1164	I = Interpolation_function(time, Q,
	1165	vertex_coordinates = points,
	1166	triangles = triangles,
	1167	interpolation_points = interpolation_points,
	1168	verbose = False)
	1169
	1170
	1171	answer = linear_function(interpolation_points)
	1172
	1173	t = time[0]
	1174	for j in range(50): #t in [1, 6]
	1175	for id in range(len(interpolation_points)):
	1176	assert allclose(I(t, id), answer[id])
	1177
	1178	t += 0.1
	1179
	1180
	1181	try:
	1182	I(1)
	1183	except:
	1184	pass
	1185	else:
	1186	raise 'Should raise exception'
	1187
	1188
	1189
[2750]	1190	def interpolation_test_interpolation_function(self):
[2526]	1191	"""Test spatio-temporal interpolation
	1192	Test that spatio temporal function performs the correct
	1193	interpolations in both time and space
	1194	"""
	1195
	1196
	1197	#Three timesteps
	1198	time = [1.0, 5.0, 6.0]
	1199
	1200
	1201	#Setup mesh used to represent fitted function
	1202	a = [0.0, 0.0]
	1203	b = [0.0, 2.0]
	1204	c = [2.0, 0.0]
	1205	d = [0.0, 4.0]
	1206	e = [2.0, 2.0]
	1207	f = [4.0, 0.0]
	1208
	1209	points = [a, b, c, d, e, f]
	1210	#bac, bce, ecf, dbe
	1211	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1212
	1213
	1214	#New datapoints where interpolated values are sought
	1215	interpolation_points = [[ 0.0, 0.0],
	1216	[ 0.5, 0.5],
	1217	[ 0.7, 0.7],
	1218	[ 1.0, 0.5],
	1219	[ 2.0, 0.4],
	1220	[ 2.8, 1.2]]
	1221
	1222
	1223	#One quantity
	1224	Q = zeros( (3,6), Float )
	1225
	1226	#Linear in time and space
	1227	for i, t in enumerate(time):
	1228	Q[i, :] = t*linear_function(points)
	1229
	1230
	1231	#Check interpolation of one quantity using interpolaton points)
	1232	I = Interpolation_function(time, Q,
	1233	vertex_coordinates = points,
	1234	triangles = triangles,
	1235	interpolation_points = interpolation_points,
	1236	verbose = False)
	1237
	1238
	1239	answer = linear_function(interpolation_points)
	1240
	1241	t = time[0]
	1242	for j in range(50): #t in [1, 6]
	1243	for id in range(len(interpolation_points)):
	1244	assert allclose(I(t, id), t*answer[id])
	1245
	1246	t += 0.1
	1247
	1248
	1249	try:
	1250	I(1)
	1251	except:
	1252	pass
	1253	else:
	1254	raise 'Should raise exception'
	1255
	1256	#
	1257	#interpolation_points = [[ 0.0, 0.0],
	1258	# [ 0.5, 0.5],
	1259	# [ 0.7, 0.7],
	1260	# [-13, 65], #Outside
	1261	# [ 1.0, 0.5],
	1262	# [ 2.0, 0.4],
	1263	# [ 2.8, 1.2]]
	1264	#
	1265	#try:
	1266	# I = Interpolation_function(time, Q,
	1267	# vertex_coordinates = points,
	1268	# triangles = triangles,
	1269	# interpolation_points = interpolation_points,
	1270	# verbose = False)
	1271	#except:
	1272	# pass
	1273	#else:
	1274	# raise 'Should raise exception'
	1275
	1276
	1277
	1278
	1279
	1280	def test_fit_and_interpolation_with_different_origins(self):
	1281	"""Fit a surface to one set of points. Then interpolate that surface
	1282	using another set of points.
	1283	This test tests situtaion where points and mesh belong to a different
	1284	coordinate system as defined by origin.
	1285	"""
	1286
	1287	#Setup mesh used to represent fitted function
	1288	a = [0.0, 0.0]
	1289	b = [0.0, 2.0]
	1290	c = [2.0, 0.0]
	1291	d = [0.0, 4.0]
	1292	e = [2.0, 2.0]
	1293	f = [4.0, 0.0]
	1294
	1295	points = [a, b, c, d, e, f]
	1296	#bac, bce, ecf, dbe, daf, dae
	1297	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1298
	1299	#Datapoints to fit from
	1300	data_points1 = [[ 0.66666667, 0.66666667],
	1301	[ 1.33333333, 1.33333333],
	1302	[ 2.66666667, 0.66666667],
	1303	[ 0.66666667, 2.66666667],
	1304	[ 0.0, 1.0],
	1305	[ 0.0, 3.0],
	1306	[ 1.0, 0.0],
	1307	[ 1.0, 1.0],
	1308	[ 1.0, 2.0],
	1309	[ 1.0, 3.0],
	1310	[ 2.0, 1.0],
	1311	[ 3.0, 0.0],
	1312	[ 3.0, 1.0]]
	1313
	1314
	1315	#First check that things are OK when using same origin
	1316	mesh_origin = (56, 290000, 618000) #zone, easting, northing
	1317	data_origin = (56, 290000, 618000) #zone, easting, northing
	1318
	1319
	1320	#Fit surface to mesh
	1321	interp = Interpolation(points, triangles, data_points1,
	1322	alpha=0.0,
	1323	data_origin = data_origin,
	1324	mesh_origin = mesh_origin)
	1325
	1326	z = linear_function(data_points1) #Example z-values
	1327	f = interp.fit(z) #Fitted values at vertices
	1328
	1329
	1330	#New datapoints where interpolated values are sought
	1331	data_points2 = [[ 0.0, 0.0],
	1332	[ 0.5, 0.5],
	1333	[ 0.7, 0.7],
	1334	[ 1.0, 0.5],
	1335	[ 2.0, 0.4],
	1336	[ 2.8, 1.2]]
	1337
	1338
	1339	#Build new A matrix based on new points
	1340	interp.build_interpolation_matrix_A(data_points2)
	1341
	1342	#Interpolate using fitted surface
[2750]	1343	#z1 = interp.interpolate(f)
[2526]	1344
	1345	#Desired result
[2750]	1346	#answer = linear_function(data_points2)
	1347	#assert allclose(z1, answer)
[2526]	1348
	1349
	1350	##############################################
	1351
	1352	#Then check situation where points are relative to a different
	1353	#origin (same zone, though, until we figure that out (FIXME))
	1354
	1355	mesh_origin = (56, 290000, 618000) #zone, easting, northing
	1356	data_origin = (56, 10000, 10000) #zone, easting, northing
	1357
	1358	#Shift datapoints according to new origin
	1359
	1360	for k in range(len(data_points1)):
	1361	data_points1[k][0] += mesh_origin[1] - data_origin[1]
	1362	data_points1[k][1] += mesh_origin[2] - data_origin[2]
	1363
	1364	for k in range(len(data_points2)):
	1365	data_points2[k][0] += mesh_origin[1] - data_origin[1]
	1366	data_points2[k][1] += mesh_origin[2] - data_origin[2]
	1367
	1368
	1369
	1370	#Fit surface to mesh
	1371	interp = Interpolation(points, triangles, data_points1,
	1372	alpha=0.0,
	1373	data_origin = data_origin,
	1374	mesh_origin = mesh_origin)
	1375
	1376	f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
	1377
	1378	assert allclose(f,f1), 'Fit should have been unaltered'
	1379
	1380
	1381	#Build new A matrix based on new points
	1382	interp.build_interpolation_matrix_A(data_points2)
	1383
	1384	#Interpolate using fitted surface
[2750]	1385	#z1 = interp.interpolate(f)
	1386	#assert allclose(z1, answer)
[2526]	1387
	1388
	1389	#########################################################
	1390	#Finally try to relate data_points2 to new origin without
	1391	#rebuilding matrix
	1392
	1393	data_origin = (56, 2000, 2000) #zone, easting, northing
	1394	for k in range(len(data_points2)):
	1395	data_points2[k][0] += 8000
	1396	data_points2[k][1] += 8000
	1397
	1398	#Build new A matrix based on new points
	1399	interp.build_interpolation_matrix_A(data_points2,
	1400	data_origin = data_origin)
	1401
	1402	#Interpolate using fitted surface
[2750]	1403	#z1 = interp.interpolate(f)
	1404	#assert allclose(z1, answer)
[2526]	1405
	1406
	1407
	1408	def test_fit_to_mesh_w_georef(self):
	1409	"""Simple check that georef works at the fit_to_mesh level
	1410	"""
	1411
	1412	from coordinate_transforms.geo_reference import Geo_reference
	1413
	1414	#Mesh
	1415	vertex_coordinates = [[0.76, 0.76],
	1416	[0.76, 5.76],
	1417	[5.76, 0.76]]
	1418	triangles = [[0,2,1]]
	1419
	1420	mesh_geo = Geo_reference(56,-0.76,-0.76)
	1421
	1422
	1423	#Data
	1424	data_points = [[ 201.0, 401.0],
	1425	[ 201.0, 403.0],
	1426	[ 203.0, 401.0]]
	1427
	1428	z = [2, 4, 4]
	1429
	1430	data_geo = Geo_reference(56,-200,-400)
	1431
	1432	#Fit
	1433	zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
	1434	data_origin = data_geo.get_origin(),
	1435	mesh_origin = mesh_geo.get_origin(),
	1436	alpha = 0)
	1437	assert allclose( zz, [0,5,5] )
	1438
	1439
	1440	def test_fit_to_mesh_file(self):
	1441	from load_mesh.loadASCII import import_mesh_file, \
	1442	export_mesh_file
	1443	import tempfile
	1444	import os
	1445
	1446	# create a .tsh file, no user outline
	1447	mesh_dic = {}
	1448	mesh_dic['vertices'] = [[0.0, 0.0],
	1449	[0.0, 5.0],
	1450	[5.0, 0.0]]
	1451	mesh_dic['triangles'] = [[0, 2, 1]]
	1452	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1453	mesh_dic['triangle_tags'] = ['']
	1454	mesh_dic['vertex_attributes'] = [[], [], []]
	1455	mesh_dic['vertiex_attribute_titles'] = []
	1456	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1457	mesh_dic['segment_tags'] = ['external',
	1458	'external',
	1459	'external']
	1460	mesh_file = tempfile.mktemp(".tsh")
	1461	export_mesh_file(mesh_file,mesh_dic)
	1462
	1463	# create an .xya file
	1464	point_file = tempfile.mktemp(".xya")
	1465	fd = open(point_file,'w')
	1466	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1467	fd.close()
	1468
	1469	mesh_output_file = tempfile.mktemp(".tsh")
	1470	fit_to_mesh_file(mesh_file,
	1471	point_file,
	1472	mesh_output_file,
	1473	alpha = 0.0)
	1474	# load in the .tsh file we just wrote
	1475	mesh_dic = import_mesh_file(mesh_output_file)
	1476	#print "mesh_dic",mesh_dic
	1477	ans =[[0.0, 0.0],
	1478	[5.0, 10.0],
	1479	[5.0,10.0]]
	1480	assert allclose(mesh_dic['vertex_attributes'],ans)
	1481
	1482	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1483	['elevation','stage'],
	1484	'test_fit_to_mesh_file failed')
	1485
	1486	#clean up
	1487	os.remove(mesh_file)
	1488	os.remove(point_file)
	1489	os.remove(mesh_output_file)
	1490
	1491	def test_fit_to_mesh_file3(self):
	1492	from load_mesh.loadASCII import import_mesh_file, \
	1493	export_mesh_file
	1494	import tempfile
	1495	import os
	1496
	1497	# create a .tsh file, no user outline
	1498	mesh_dic = {}
	1499	mesh_dic['vertices'] = [[0.76, 0.76],
	1500	[0.76, 5.76],
	1501	[5.76, 0.76]]
	1502	mesh_dic['triangles'] = [[0, 2, 1]]
	1503	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1504	mesh_dic['triangle_tags'] = ['']
	1505	mesh_dic['vertex_attributes'] = [[], [], []]
	1506	mesh_dic['vertiex_attribute_titles'] = []
	1507	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1508	mesh_dic['segment_tags'] = ['external',
	1509	'external',
	1510	'external']
	1511	mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
	1512	mesh_file = tempfile.mktemp(".tsh")
	1513	export_mesh_file(mesh_file,mesh_dic)
	1514
	1515	# create an .xya file
	1516	point_file = tempfile.mktemp(".xya")
	1517	fd = open(point_file,'w')
	1518	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1519	fd.close()
	1520
	1521	mesh_output_file = tempfile.mktemp(".tsh")
	1522	fit_to_mesh_file(mesh_file,
	1523	point_file,
	1524	mesh_output_file,
	1525	alpha = 0.0)
	1526	# load in the .tsh file we just wrote
	1527	mesh_dic = import_mesh_file(mesh_output_file)
	1528	#print "mesh_dic",mesh_dic
	1529	ans =[[0.0, 0.0],
	1530	[5.0, 10.0],
	1531	[5.0,10.0]]
	1532	assert allclose(mesh_dic['vertex_attributes'],ans)
	1533
	1534	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1535	['elevation','stage'],
	1536	'test_fit_to_mesh_file failed')
	1537
	1538	#clean up
	1539	os.remove(mesh_file)
	1540	os.remove(point_file)
	1541	os.remove(mesh_output_file)
	1542
	1543	def test_fit_to_mesh_file4(self):
	1544	from load_mesh.loadASCII import import_mesh_file, \
	1545	export_mesh_file
	1546	import tempfile
	1547	import os
	1548
	1549	# create a .tsh file, no user outline
	1550	mesh_dic = {}
	1551	mesh_dic['vertices'] = [[0.76, 0.76],
	1552	[0.76, 5.76],
	1553	[5.76, 0.76]]
	1554	mesh_dic['triangles'] = [[0, 2, 1]]
	1555	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1556	mesh_dic['triangle_tags'] = ['']
	1557	mesh_dic['vertex_attributes'] = [[], [], []]
	1558	mesh_dic['vertiex_attribute_titles'] = []
	1559	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1560	mesh_dic['segment_tags'] = ['external',
	1561	'external',
	1562	'external']
	1563	mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
	1564	mesh_file = tempfile.mktemp(".tsh")
	1565	export_mesh_file(mesh_file,mesh_dic)
	1566
	1567	geo_ref = Geo_reference(56,-200,-400)
	1568	# create an .xya file
	1569	point_file = tempfile.mktemp(".xya")
	1570	fd = open(point_file,'w')
	1571	fd.write("elevation, stage \n 201.0, 401.0,2.,4 \n 201.0, 403.0,4,8 \n 203.0, 401.0,4.,8 \n")
	1572	geo_ref.write_ASCII(fd)
	1573	fd.close()
	1574
	1575	mesh_output_file = tempfile.mktemp(".tsh")
	1576	fit_to_mesh_file(mesh_file,
	1577	point_file,
	1578	mesh_output_file,
	1579	alpha = 0.0)
	1580	# load in the .tsh file we just wrote
	1581	mesh_dic = import_mesh_file(mesh_output_file)
	1582	#print "mesh_dic",mesh_dic
	1583	ans =[[0.0, 0.0],
	1584	[5.0, 10.0],
	1585	[5.0, 10.0]]
	1586	assert allclose(mesh_dic['vertex_attributes'],ans)
	1587
	1588	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1589	['elevation','stage'],
	1590	'test_fit_to_mesh_file failed')
	1591
	1592	#clean up
	1593	os.remove(mesh_file)
	1594	os.remove(point_file)
	1595	os.remove(mesh_output_file)
	1596
	1597	def test_fit_to_mesh_fileII(self):
	1598	from load_mesh.loadASCII import import_mesh_file, \
	1599	export_mesh_file
	1600	import tempfile
	1601	import os
	1602
	1603	# create a .tsh file, no user outline
	1604	mesh_dic = {}
	1605	mesh_dic['vertices'] = [[0.0, 0.0],
	1606	[0.0, 5.0],
	1607	[5.0, 0.0]]
	1608	mesh_dic['triangles'] = [[0, 2, 1]]
	1609	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1610	mesh_dic['triangle_tags'] = ['']
	1611	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1612	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1613	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1614	mesh_dic['segment_tags'] = ['external',
	1615	'external',
	1616	'external']
	1617	mesh_file = tempfile.mktemp(".tsh")
	1618	export_mesh_file(mesh_file,mesh_dic)
	1619
	1620	# create an .xya file
	1621	point_file = tempfile.mktemp(".xya")
	1622	fd = open(point_file,'w')
	1623	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1624	fd.close()
	1625
	1626	mesh_output_file = "new_triangle.tsh"
	1627	fit_to_mesh_file(mesh_file,
	1628	point_file,
	1629	mesh_output_file,
	1630	alpha = 0.0)
	1631	# load in the .tsh file we just wrote
	1632	mesh_dic = import_mesh_file(mesh_output_file)
	1633
	1634	assert allclose(mesh_dic['vertex_attributes'],
	1635	[[1.0, 2.0,0.0, 0.0],
	1636	[1.0, 2.0,5.0, 10.0],
	1637	[1.0, 2.0,5.0,10.0]])
	1638
	1639	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1640	['density', 'temp','elevation','stage'],
	1641	'test_fit_to_mesh_file failed')
	1642
	1643	#clean up
	1644	os.remove(mesh_file)
	1645	os.remove(mesh_output_file)
	1646	os.remove(point_file)
	1647
	1648	def test_fit_to_mesh_file_errors(self):
	1649	from load_mesh.loadASCII import import_mesh_file, export_mesh_file
	1650	import tempfile
	1651	import os
	1652
	1653	# create a .tsh file, no user outline
	1654	mesh_dic = {}
	1655	mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
	1656	mesh_dic['triangles'] = [[0, 2, 1]]
	1657	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1658	mesh_dic['triangle_tags'] = ['']
	1659	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1660	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1661	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1662	mesh_dic['segment_tags'] = ['external', 'external','external']
	1663	mesh_file = tempfile.mktemp(".tsh")
	1664	export_mesh_file(mesh_file,mesh_dic)
	1665
	1666	# create an .xya file
	1667	point_file = tempfile.mktemp(".xya")
	1668	fd = open(point_file,'w')
	1669	fd.write("elevation stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1670	fd.close()
	1671
	1672	mesh_output_file = "new_triangle.tsh"
	1673	try:
	1674	fit_to_mesh_file(mesh_file, point_file,
	1675	mesh_output_file, display_errors = False)
	1676	except IOError:
	1677	pass
	1678	else:
	1679	#self.failUnless(0 ==1, 'Bad file did not raise error!')
	1680	raise 'Bad file did not raise error!'
	1681
	1682	#clean up
	1683	os.remove(mesh_file)
	1684	os.remove(point_file)
	1685
	1686	def test_fit_to_mesh_file_errorsII(self):
	1687	from load_mesh.loadASCII import import_mesh_file, export_mesh_file
	1688	import tempfile
	1689	import os
	1690
	1691	# create a .tsh file, no user outline
	1692	mesh_file = tempfile.mktemp(".tsh")
	1693	fd = open(mesh_file,'w')
	1694	fd.write("unit testing a bad .tsh file \n")
	1695	fd.close()
	1696
	1697	# create an .xya file
	1698	point_file = tempfile.mktemp(".xya")
	1699	fd = open(point_file,'w')
	1700	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1701	fd.close()
	1702
	1703	mesh_output_file = "new_triangle.tsh"
	1704	try:
	1705	fit_to_mesh_file(mesh_file, point_file,
	1706	mesh_output_file, display_errors = False)
	1707	except IOError:
	1708	pass
	1709	else:
	1710	raise 'Bad file did not raise error!'
	1711
	1712	#clean up
	1713	os.remove(mesh_file)
	1714	os.remove(point_file)
	1715
	1716	def test_fit_to_mesh_file_errorsIII(self):
	1717	from load_mesh.loadASCII import import_mesh_file, export_mesh_file
	1718	import tempfile
	1719	import os
	1720
	1721	# create a .tsh file, no user outline
	1722	mesh_dic = {}
	1723	mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
	1724	mesh_dic['triangles'] = [[0, 2, 1]]
	1725	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1726	mesh_dic['triangle_tags'] = ['']
	1727	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1728	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1729	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1730	mesh_dic['segment_tags'] = ['external', 'external','external']
	1731	mesh_file = tempfile.mktemp(".tsh")
	1732	export_mesh_file(mesh_file,mesh_dic)
	1733
	1734	# create an .xya file
	1735	point_file = tempfile.mktemp(".xya")
	1736	fd = open(point_file,'w')
	1737	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1738	fd.close()
	1739
	1740	#This a deliberately illegal filename to invoke the error.
	1741	mesh_output_file = ".../\z\z:ya.tsh"
	1742
	1743	try:
	1744	fit_to_mesh_file(mesh_file, point_file,
	1745	mesh_output_file, display_errors = False)
	1746	except IOError:
	1747	pass
	1748	else:
	1749	raise 'Bad file did not raise error!'
	1750
	1751	#clean up
	1752	os.remove(mesh_file)
	1753	os.remove(point_file)
	1754
	1755	## FIXME? Running from the Comand line isn't in vogue these days
	1756	# The test was breaking when test_all at the inundation level was running
	1757	# was running it.issue - not running the test in this directory
	1758	def Bad_test_fit_to_mesh_file_errorsIV(self):
	1759	import os
	1760	command = '%s least_squares.py q q q e n 0.9 n' %(sys.executable)
	1761	status = os.system(command)
	1762	self.failUnless(status%255 == 1,
	1763	'command prompt least_squares.py failed. Incorect exit status.')
	1764
	1765	def test_fit_to_msh_netcdf_fileII(self):
	1766	from load_mesh.loadASCII import import_mesh_file, export_mesh_file
	1767	import tempfile
	1768	import os
	1769
	1770	# create a .tsh file, no user outline
	1771	mesh_dic = {}
	1772	mesh_dic['vertices'] = [[0.0, 0.0],
	1773	[0.0, 5.0],
	1774	[5.0, 0.0]]
	1775	mesh_dic['triangles'] = [[0, 2, 1]]
	1776	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1777	mesh_dic['triangle_tags'] = ['']
	1778	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1779	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1780	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1781	mesh_dic['segment_tags'] = ['external',
	1782	'external',
	1783	'external']
	1784	mesh_file = tempfile.mktemp(".msh")
	1785	export_mesh_file(mesh_file,mesh_dic)
	1786
	1787	# create an .xya file
	1788	point_file = tempfile.mktemp(".xya")
	1789	fd = open(point_file,'w')
	1790	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1791	fd.close()
	1792
	1793	mesh_output_file = "new_triangle.msh"
	1794	fit_to_mesh_file(mesh_file,
	1795	point_file,
	1796	mesh_output_file,
	1797	alpha = 0.0)
	1798	# load in the .tsh file we just wrote
	1799	mesh_dic = import_mesh_file(mesh_output_file)
	1800
	1801	assert allclose(mesh_dic['vertex_attributes'],
	1802	[[1.0, 2.0,0.0, 0.0],
	1803	[1.0, 2.0,5.0, 10.0],
	1804	[1.0, 2.0,5.0,10.0]])
	1805
	1806	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1807	['density', 'temp','elevation','stage'],
	1808	'test_fit_to_mesh_file failed')
	1809
	1810	#clean up
	1811	os.remove(mesh_file)
	1812	os.remove(mesh_output_file)
	1813	os.remove(point_file)
	1814
[2585]	1815
	1816
	1817	def test_fit_using_fit_to_mesh(self):
	1818	"""Fit a surface to one set of points. Then interpolate that surface
	1819	using another set of points.
	1820	"""
	1821
	1822	#Setup mesh used to represent fitted function
	1823	a = [0.0, 0.0]
	1824	b = [0.0, 2.0]
	1825	c = [2.0, 0.0]
	1826	d = [0.0, 4.0]
	1827	e = [2.0, 2.0]
	1828	f = [4.0, 0.0]
	1829
	1830	points = [a, b, c, d, e, f]
	1831	#bac, bce, ecf, dbe, daf, dae
	1832	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1833
	1834	#Datapoints to fit from
	1835	data_points1 = [[ 0.66666667, 0.66666667],
	1836	[ 1.33333333, 1.33333333],
	1837	[ 2.66666667, 0.66666667],
	1838	[ 0.66666667, 2.66666667],
	1839	[ 0.0, 1.0],
	1840	[ 0.0, 3.0],
	1841	[ 1.0, 0.0],
	1842	[ 1.0, 1.0],
	1843	[ 15, -17], #Outside mesh
	1844	[ 1.0, 2.0],
	1845	[ 1.0, 3.0],
	1846	[ 2.0, 1.0],
	1847	[ 3.0, 0.0],
	1848	[ 3.0, 1.0]]
	1849
	1850	#Fit surface to mesh
	1851	z = linear_function(data_points1) #Example z-values
[2666]	1852	v = fit_to_mesh(points, triangles, data_points1, z, alpha=0.0,
	1853	precrop=True, verbose=False)
	1854
[2585]	1855	assert allclose(linear_function(points), v)
	1856
[2666]	1857
[2585]	1858
[2666]	1859	def test_acceptable_overshoot(self):
	1860	"""Fit a surface to one set of points. Then interpolate that surface
	1861	using another set of points.
	1862	Check that exceedance in fitted values are caught.
	1863	"""
	1864
	1865	#Setup mesh used to represent fitted function
	1866	a = [0.0, 0.0]
	1867	b = [0.0, 2.0]
	1868	c = [2.0, 0.0]
	1869	d = [0.0, 4.0]
	1870	e = [2.0, 2.0]
	1871	f = [4.0, 0.0]
	1872
	1873	points = [a, b, c, d, e, f]
	1874	#bac, bce, ecf, dbe, daf, dae
	1875	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1876
	1877	#Datapoints to fit from
	1878	data_points1 = [[ 0.66666667, 0.66666667],
	1879	[ 1.33333333, 1.33333333],
	1880	[ 2.66666667, 0.66666667],
	1881	[ 0.66666667, 2.66666667],
	1882	[ 0.0, 1.0],
	1883	[ 0.0, 3.0],
	1884	[ 1.0, 0.0],
	1885	[ 1.0, 1.0],
	1886	[ 15, -17], #Outside mesh
	1887	[ 1.0, 2.0],
	1888	[ 1.0, 3.0],
	1889	[ 2.0, 1.0],
	1890	[ 3.0, 0.0],
	1891	[ 3.0, 1.0]]
	1892
	1893	#Fit surface to mesh
	1894	z = linear_function(data_points1) #Example z-values
	1895
	1896	try:
	1897	v = fit_to_mesh(points, triangles, data_points1, z, alpha=0.0,
	1898	acceptable_overshoot = 0.2,
	1899	precrop=True, verbose=False)
	1900	except FittingError, e:
	1901	pass
	1902	else:
	1903	raise 'Should have raised exception'
	1904
	1905
	1906	#assert allclose(linear_function(points), v)
	1907
	1908
	1909
[2526]	1910	#-------------------------------------------------------------
	1911	if __name__ == "__main__":
[2665]	1912	#suite = unittest.makeSuite(Test_Least_Squares,'test_smooth_attributes_to_mesh_function')
[2778]	1913	#suite = unittest.makeSuite(Test_Least_Squares,'test_datapoint_at_centroid')
[2666]	1914	suite = unittest.makeSuite(Test_Least_Squares,'test')
[2526]	1915
	1916	#suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
	1917	#suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
	1918	runner = unittest.TextTestRunner(verbosity=1)
	1919	runner.run(suite)
	1920
	1921
	1922
	1923
	1924

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

Context Navigation

source: inundation/pyvolution/test_least_squares.py @ 3236

Download in other formats: