Context Navigation

source: inundation/fit_interpolate/interpolate.py @ 3409

Last change on this file since 3409 was 3409, checked in by duncan, 18 years ago
comments
File size: 26.9 KB

Line
1	"""Least squares interpolation.
2
3	Implements a least-squares interpolation.
4	Putting mesh data onto points.
5
6	Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
7	Geoscience Australia, 2004.
8
9	DESIGN ISSUES
10	* what variables should be global?
11	- if there are no global vars functions can be moved around alot easier
12
13	* The public interface to Interpolate
14	__init__
15	interpolate
16	interpolate_block
17
18	"""
19
20	import time
21	import os
22	from warnings import warn
23
24	from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
25	ArrayType, allclose, take, NewAxis, arange
26
27	from caching.caching import cache
28	from pyvolution.neighbour_mesh import Mesh
29	from utilities.sparse import Sparse, Sparse_CSR
30	from utilities.cg_solve import conjugate_gradient, VectorShapeError
31	from coordinate_transforms.geo_reference import Geo_reference
32	from pyvolution.quad import build_quadtree
33	from utilities.numerical_tools import ensure_numeric, mean, INF
34	from utilities.polygon import in_and_outside_polygon
35	from geospatial_data.geospatial_data import Geospatial_data
36	from fit_interpolate.search_functions import search_tree_of_vertices
37	from fit_interpolate.general_fit_interpolate import FitInterpolate
38
39
40
41	class Interpolate (FitInterpolate):
42
43	def __init__(self,
44	vertex_coordinates,
45	triangles,
46	mesh_origin=None,
47	verbose=False,
48	max_vertices_per_cell=30):
49
50
51	""" Build interpolation matrix mapping from
52	function values at vertices to function values at data points
53
54	Inputs:
55
56	vertex_coordinates: List of coordinate pairs [xi, eta] of
57	points constituting a mesh (or an m x 2 Numeric array or
58	a geospatial object)
59	Points may appear multiple times
60	(e.g. if vertices have discontinuities)
61
62	triangles: List of 3-tuples (or a Numeric array) of
63	integers representing indices of all vertices in the mesh.
64
65	mesh_origin: A geo_reference object or 3-tuples consisting of
66	UTM zone, easting and northing.
67	If specified vertex coordinates are assumed to be
68	relative to their respective origins.
69
70	max_vertices_per_cell: Number of vertices in a quad tree cell
71	at which the cell is split into 4.
72
73	Note: Don't supply a vertex coords as a geospatial object and
74	a mesh origin, since geospatial has its own mesh origin.
75	"""
76
77	# Initialise variabels
78	self._A_can_be_reused = False
79	self._point_coordinates = None
80
81	FitInterpolate.__init__(self,
82	vertex_coordinates,
83	triangles,
84	mesh_origin,
85	verbose,
86	max_vertices_per_cell)
87
88	# FIXME: What is a good start_blocking_len value?
89	def interpolate(self, f, point_coordinates = None,
90	start_blocking_len = 500000, verbose=False):
91	"""Interpolate mesh data f to determine values, z, at points.
92
93	f is the data on the mesh vertices.
94
95	The mesh values representing a smooth surface are
96	assumed to be specified in f.
97
98	Inputs:
99	f: Vector or array of data at the mesh vertices.
100	If f is an array, interpolation will be done for each column as
101	per underlying matrix-matrix multiplication
102	point_coordinates: Interpolate mesh data to these positions.
103	List of coordinate pairs [x, y] of
104	data points or an nx2 Numeric array or a Geospatial_data object
105
106	If point_coordinates is absent, the points inputted last time
107	this method was called are used, if possible.
108	start_blocking_len: If the # of points is more or greater than this,
109	start blocking
110
111	Output:
112	Interpolated values at inputted points (z).
113	"""
114	#print "point_coordinates interpolate.interpolate",point_coordinates
115	if isinstance(point_coordinates,Geospatial_data):
116	point_coordinates = point_coordinates.get_data_points( \
117	absolute = True)
118
119	# Can I interpolate, based on previous point_coordinates?
120	if point_coordinates is None:
121	if self._A_can_be_reused is True and \
122	len(self._point_coordinates) < start_blocking_len:
123	z = self._get_point_data_z(f,
124	verbose=verbose)
125	elif self._point_coordinates is not None:
126	# if verbose, give warning
127	if verbose:
128	print 'WARNING: Recalculating A matrix, due to blocking.'
129	point_coordinates = self._point_coordinates
130	else:
131	#There are no good point_coordinates. import sys; sys.exit()
132	msg = 'ERROR (interpolate.py): No point_coordinates inputted'
133	raise Exception(msg)
134
135
136	if point_coordinates is not None:
137	self._point_coordinates = point_coordinates
138	if len(point_coordinates) < start_blocking_len or \
139	start_blocking_len == 0:
140	self._A_can_be_reused = True
141	z = self.interpolate_block(f, point_coordinates,
142	verbose=verbose)
143	else:
144	#Handle blocking
145	self._A_can_be_reused = False
146	start=0
147	# creating a dummy array to concatenate to.
148
149	f = ensure_numeric(f, Float)
150	#print "f.shape",f.shape
151	if len(f.shape) > 1:
152	z = zeros((0,f.shape[1]))
153	else:
154	z = zeros((0,))
155
156	for end in range(start_blocking_len
157	,len(point_coordinates)
158	,start_blocking_len):
159	t = self.interpolate_block(f, point_coordinates[start:end],
160	verbose=verbose)
161	#print "t", t
162	#print "z", z
163	z = concatenate((z,t))
164	start = end
165	end = len(point_coordinates)
166	t = self.interpolate_block(f, point_coordinates[start:end],
167	verbose=verbose)
168	z = concatenate((z,t))
169	return z
170
171	def interpolate_block(self, f, point_coordinates = None, verbose=False):
172	"""
173	Call this if you want to control the blocking or make sure blocking
174	doesn't occur.
175
176	Return the point data, z.
177
178	See interpolate for doc info.
179	"""
180	if isinstance(point_coordinates,Geospatial_data):
181	point_coordinates = point_coordinates.get_data_points( \
182	absolute = True)
183	if point_coordinates is not None:
184	self._A =self._build_interpolation_matrix_A(point_coordinates,
185	verbose=verbose)
186	return self._get_point_data_z(f)
187
188	def _get_point_data_z(self, f, verbose=False):
189	"""
190	Return the point data, z.
191
192	Precondition,
193	The _A matrix has been created
194	"""
195	z = self._A * f
196	# Taking into account points outside the mesh.
197	#print "self.outside_poly_indices", self.outside_poly_indices
198	#print "self.inside_poly_indices", self.inside_poly_indices
199	#print "z", z
200	for i in self.outside_poly_indices:
201	z[i] = INF
202	return z
203
204
205	def _build_interpolation_matrix_A(self,
206	point_coordinates,
207	verbose = False):
208	"""Build n x m interpolation matrix, where
209	n is the number of data points and
210	m is the number of basis functions phi_k (one per vertex)
211
212	This algorithm uses a quad tree data structure for fast binning
213	of data points
214	origin is a 3-tuple consisting of UTM zone, easting and northing.
215	If specified coordinates are assumed to be relative to this origin.
216
217	This one will override any data_origin that may be specified in
218	instance interpolation
219
220	Preconditions
221	Point_coordindates and mesh vertices have the same origin.
222	"""
223
224
225
226	#Convert point_coordinates to Numeric arrays, in case it was a list.
227	point_coordinates = ensure_numeric(point_coordinates, Float)
228
229	if verbose: print 'Getting indices inside mesh boundary'
230	#print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon()
231	self.inside_poly_indices, self.outside_poly_indices = \
232	in_and_outside_polygon(point_coordinates,
233	self.mesh.get_boundary_polygon(),
234	closed = True, verbose = verbose)
235	#print "self.inside_poly_indices",self.inside_poly_indices
236	#print "self.outside_poly_indices",self.outside_poly_indices
237	#Build n x m interpolation matrix
238	if verbose and len(self.outside_poly_indices) >0:
239	print '\n WARNING: Points outside mesh boundary. \n'
240	# Since you can block, throw a warning, not an error.
241	if verbose and 0 == len(self.inside_poly_indices):
242	print '\n WARNING: No points within the mesh! \n'
243
244	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
245	n = point_coordinates.shape[0] #Nbr of data points
246
247	if verbose: print 'Number of datapoints: %d' %n
248	if verbose: print 'Number of basis functions: %d' %m
249
250	A = Sparse(n,m)
251
252	n = len(self.inside_poly_indices)
253	#Compute matrix elements for points inside the mesh
254	for i in self.inside_poly_indices:
255	#For each data_coordinate point
256	if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
257	x = point_coordinates[i]
258	element_found, sigma0, sigma1, sigma2, k = \
259	search_tree_of_vertices(self.root, self.mesh, x)
260	#Update interpolation matrix A if necessary
261	if element_found is True:
262	#Assign values to matrix A
263
264	j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
265	j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
266	j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
267
268	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
269	js = [j0,j1,j2]
270
271	for j in js:
272	A[i,j] = sigmas[j]
273	else:
274	msg = 'Could not find triangle for point', x
275	raise Exception(msg)
276	return A
277
278	class Interpolation_function:
279	"""Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
280	which is interpolated from time series defined at vertices of
281	triangular mesh (such as those stored in sww files)
282
283	Let m be the number of vertices, n the number of triangles
284	and p the number of timesteps.
285
286	Mandatory input
287	time: px1 array of monotonously increasing times (Float)
288	quantities: Dictionary of arrays or 1 array (Float)
289	The arrays must either have dimensions pxm or mx1.
290	The resulting function will be time dependent in
291	the former case while it will be constan with
292	respect to time in the latter case.
293
294	Optional input:
295	quantity_names: List of keys into the quantities dictionary
296	vertex_coordinates: mx2 array of coordinates (Float)
297	triangles: nx3 array of indices into vertex_coordinates (Int)
298	interpolation_points: Nx2 array of coordinates to be interpolated to
299	verbose: Level of reporting
300
301
302	The quantities returned by the callable object are specified by
303	the list quantities which must contain the names of the
304	quantities to be returned and also reflect the order, e.g. for
305	the shallow water wave equation, on would have
306	quantities = ['stage', 'xmomentum', 'ymomentum']
307
308	The parameter interpolation_points decides at which points interpolated
309	quantities are to be computed whenever object is called.
310	If None, return average value
311	"""
312
313
314
315	def __init__(self,
316	time,
317	quantities,
318	quantity_names = None,
319	vertex_coordinates = None,
320	triangles = None,
321	interpolation_points = None,
322	verbose = False):
323	"""Initialise object and build spatial interpolation if required
324	"""
325
326	from Numeric import array, zeros, Float, alltrue, concatenate,\
327	reshape, ArrayType
328
329
330	#from util import mean, ensure_numeric
331	from config import time_format
332	import types
333
334
335	#Check temporal info
336	time = ensure_numeric(time)
337	msg = 'Time must be a monotonuosly '
338	msg += 'increasing sequence %s' %time
339	assert alltrue(time[1:] - time[:-1] >= 0 ), msg
340
341
342	#Check if quantities is a single array only
343	if type(quantities) != types.DictType:
344	quantities = ensure_numeric(quantities)
345	quantity_names = ['Attribute']
346
347	#Make it a dictionary
348	quantities = {quantity_names[0]: quantities}
349
350
351	#Use keys if no names are specified
352	if quantity_names is None:
353	quantity_names = quantities.keys()
354
355
356	#Check spatial info
357	if vertex_coordinates is None:
358	self.spatial = False
359	else:
360	vertex_coordinates = ensure_numeric(vertex_coordinates)
361
362	assert triangles is not None, 'Triangles array must be specified'
363	triangles = ensure_numeric(triangles)
364	self.spatial = True
365
366
367	#Save for use with statistics
368
369	self.quantities_range = {}
370	for name in quantity_names:
371	q = quantities[name][:].flat
372	self.quantities_range[name] = [min(q), max(q)]
373
374	self.quantity_names = quantity_names
375	#self.quantities = quantities #Takes too much memory
376	self.vertex_coordinates = vertex_coordinates
377	self.interpolation_points = interpolation_points
378	self.time = time[:] # Time assumed to be relative to starttime
379	self.index = 0 # Initial time index
380	self.precomputed_values = {}
381
382
383
384
385
386
387	#Precomputed spatial interpolation if requested
388	if interpolation_points is not None:
389	if self.spatial is False:
390	raise 'Triangles and vertex_coordinates must be specified'
391
392	try:
393	self.interpolation_points = ensure_numeric(interpolation_points)
394	except:
395	msg = 'Interpolation points must be an N x 2 Numeric array '+\
396	'or a list of points\n'
397	msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
398	'...')
399	raise msg
400
401
402	m = len(self.interpolation_points)
403	p = len(self.time)
404
405	for name in quantity_names:
406	self.precomputed_values[name] = zeros((p, m), Float)
407
408	#Build interpolator
409	interpol = Interpolate(vertex_coordinates,
410	triangles,
411	#point_coordinates = \
412	#self.interpolation_points,
413	#alpha = 0,
414	verbose = verbose)
415
416	if verbose: print 'Interpolate'
417	for i, t in enumerate(self.time):
418	#Interpolate quantities at this timestep
419	if verbose and i%((p+10)/10)==0:
420	print ' time step %d of %d' %(i, p)
421
422	for name in quantity_names:
423	if len(quantities[name].shape) == 2:
424	result = interpol.interpolate(quantities[name][i,:],
425	point_coordinates = \
426	self.interpolation_points)
427	else:
428	#Assume no time dependency
429	result = interpol.interpolate(quantities[name][:],
430	point_coordinates = \
431	self.interpolation_points)
432
433	self.precomputed_values[name][i, :] = result
434
435	#Report
436	if verbose:
437	print self.statistics()
438	#self.print_statistics()
439
440	else:
441	#Store quantitites as is
442	for name in quantity_names:
443	self.precomputed_values[name] = quantities[name]
444
445
446	def __repr__(self):
447	#return 'Interpolation function (spatio-temporal)'
448	return self.statistics()
449
450
451	def __call__(self, t, point_id = None, x = None, y = None):
452	"""Evaluate f(t), f(t, point_id) or f(t, x, y)
453
454	Inputs:
455	t: time - Model time. Must lie within existing timesteps
456	point_id: index of one of the preprocessed points.
457	x, y: Overrides location, point_id ignored
458
459	If spatial info is present and all of x,y,point_id
460	are None an exception is raised
461
462	If no spatial info is present, point_id and x,y arguments are ignored
463	making f a function of time only.
464
465
466	FIXME: point_id could also be a slice
467	FIXME: What if x and y are vectors?
468	FIXME: What about f(x,y) without t?
469	"""
470
471	from math import pi, cos, sin, sqrt
472	from Numeric import zeros, Float
473	from util import mean
474
475	if self.spatial is True:
476	if point_id is None:
477	if x is None or y is None:
478	msg = 'Either point_id or x and y must be specified'
479	raise Exception(msg)
480	else:
481	if self.interpolation_points is None:
482	msg = 'Interpolation_function must be instantiated ' +\
483	'with a list of interpolation points before parameter ' +\
484	'point_id can be used'
485	raise Exception(msg)
486
487	msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
488	msg += ' does not match model time: %.16f\n' %t
489	if t < self.time[0]: raise Exception(msg)
490	if t > self.time[-1]: raise Exception(msg)
491
492	oldindex = self.index #Time index
493
494	#Find current time slot
495	while t > self.time[self.index]: self.index += 1
496	while t < self.time[self.index]: self.index -= 1
497
498	if t == self.time[self.index]:
499	#Protect against case where t == T[-1] (last time)
500	# - also works in general when t == T[i]
501	ratio = 0
502	else:
503	#t is now between index and index+1
504	ratio = (t - self.time[self.index])/\
505	(self.time[self.index+1] - self.time[self.index])
506
507	#Compute interpolated values
508	q = zeros(len(self.quantity_names), Float)
509	#print "self.precomputed_values", self.precomputed_values
510	for i, name in enumerate(self.quantity_names):
511	Q = self.precomputed_values[name]
512
513	if self.spatial is False:
514	#If there is no spatial info
515	assert len(Q.shape) == 1
516
517	Q0 = Q[self.index]
518	if ratio > 0: Q1 = Q[self.index+1]
519
520	else:
521	if x is not None and y is not None:
522	#Interpolate to x, y
523
524	raise 'x,y interpolation not yet implemented'
525	else:
526	#Use precomputed point
527	Q0 = Q[self.index, point_id]
528	if ratio > 0:
529	Q1 = Q[self.index+1, point_id]
530
531	#Linear temporal interpolation
532	if ratio > 0:
533	if Q0 == INF and Q1 == INF:
534	q[i] = Q0
535	else:
536	q[i] = Q0 + ratio*(Q1 - Q0)
537	else:
538	q[i] = Q0
539
540
541	#Return vector of interpolated values
542	#if len(q) == 1:
543	# return q[0]
544	#else:
545	# return q
546
547
548	#Return vector of interpolated values
549	#FIXME:
550	if self.spatial is True:
551	return q
552	else:
553	#Replicate q according to x and y
554	#This is e.g used for Wind_stress
555	if x is None or y is None:
556	return q
557	else:
558	try:
559	N = len(x)
560	except:
561	return q
562	else:
563	from Numeric import ones, Float
564	#x is a vector - Create one constant column for each value
565	N = len(x)
566	assert len(y) == N, 'x and y must have same length'
567	res = []
568	for col in q:
569	res.append(col*ones(N, Float))
570
571	return res
572
573
574	def get_time(self):
575	"""Return model time as a vector of timesteps
576	"""
577	return self.time
578
579	def statistics(self):
580	"""Output statistics about interpolation_function
581	"""
582
583	vertex_coordinates = self.vertex_coordinates
584	interpolation_points = self.interpolation_points
585	quantity_names = self.quantity_names
586	#quantities = self.quantities
587	precomputed_values = self.precomputed_values
588
589	x = vertex_coordinates[:,0]
590	y = vertex_coordinates[:,1]
591
592	str = '------------------------------------------------\n'
593	str += 'Interpolation_function (spatio-temporal) statistics:\n'
594	str += ' Extent:\n'
595	str += ' x in [%f, %f], len(x) == %d\n'\
596	%(min(x), max(x), len(x))
597	str += ' y in [%f, %f], len(y) == %d\n'\
598	%(min(y), max(y), len(y))
599	str += ' t in [%f, %f], len(t) == %d\n'\
600	%(min(self.time), max(self.time), len(self.time))
601	str += ' Quantities:\n'
602	for name in quantity_names:
603	minq, maxq = self.quantities_range[name]
604	str += ' %s in [%f, %f]\n' %(name, minq, maxq)
605	#q = quantities[name][:].flat
606	#str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
607
608	if interpolation_points is not None:
609	str += ' Interpolation points (xi, eta):'\
610	' number of points == %d\n' %interpolation_points.shape[0]
611	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
612	max(interpolation_points[:,0]))
613	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
614	max(interpolation_points[:,1]))
615	str += ' Interpolated quantities (over all timesteps):\n'
616
617	for name in quantity_names:
618	q = precomputed_values[name][:].flat
619	str += ' %s at interpolation points in [%f, %f]\n'\
620	%(name, min(q), max(q))
621	str += '------------------------------------------------\n'
622
623	return str
624
625	def interpolate_sww(sww_file, time, interpolation_points,
626	quantity_names = None, verbose = False):
627	"""
628	obsolete.
629	use file_function in utils
630	"""
631	#open sww file
632	x, y, volumes, time, quantities = read_sww(sww_file)
633	print "x",x
634	print "y",y
635
636	print "time", time
637	print "quantities", quantities
638
639	#Add the x and y together
640	vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
641
642	#Will return the quantity values at the specified times and locations
643	interp = Interpolation_interface(
644	time,
645	quantities,
646	quantity_names = quantity_names,
647	vertex_coordinates = vertex_coordinates,
648	triangles = volumes,
649	interpolation_points = interpolation_points,
650	verbose = verbose)
651
652
653	def read_sww(file_name):
654	"""
655	obsolete - Nothing should be calling this
656
657	Read in an sww file.
658
659	Input;
660	file_name - the sww file
661
662	Output;
663	x - Vector of x values
664	y - Vector of y values
665	z - Vector of bed elevation
666	volumes - Array. Each row has 3 values, representing
667	the vertices that define the volume
668	time - Vector of the times where there is stage information
669	stage - array with respect to time and vertices (x,y)
670	"""
671
672	#FIXME Have this reader as part of data_manager?
673
674	from Scientific.IO.NetCDF import NetCDFFile
675	import tempfile
676	import sys
677	import os
678
679	#Check contents
680	#Get NetCDF
681
682	# see if the file is there. Throw a QUIET IO error if it isn't
683	# I don't think I could implement the above
684
685	#throws prints to screen if file not present
686	junk = tempfile.mktemp(".txt")
687	fd = open(junk,'w')
688	stdout = sys.stdout
689	sys.stdout = fd
690	fid = NetCDFFile(file_name, 'r')
691	sys.stdout = stdout
692	fd.close()
693	#clean up
694	os.remove(junk)
695
696	# Get the variables
697	x = fid.variables['x'][:]
698	y = fid.variables['y'][:]
699	volumes = fid.variables['volumes'][:]
700	time = fid.variables['time'][:]
701
702	keys = fid.variables.keys()
703	keys.remove("x")
704	keys.remove("y")
705	keys.remove("volumes")
706	keys.remove("time")
707	#Turn NetCDF objects into Numeric arrays
708	quantities = {}
709	for name in keys:
710	quantities[name] = fid.variables[name][:]
711
712	fid.close()
713	return x, y, volumes, time, quantities
714
715	#-------------------------------------------------------------
716	if __name__ == "__main__":
717	a = [0.0, 0.0]
718	b = [0.0, 2.0]
719	c = [2.0,0.0]
720	points = [a, b, c]
721	vertices = [ [1,0,2] ] #bac
722
723	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
724
725	interp = Interpolate(points, vertices) #, data)
726	A = interp._build_interpolation_matrix_A(data, verbose=True)
727	A = A.todense()
728	print "A",A
729	assert allclose(A, [[1./3, 1./3, 1./3]])
730	print "finished"

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

Download in other formats: