Context Navigation

source: inundation/fit_interpolate/interpolate.py @ 3330

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

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

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

Download in other formats: