Context Navigation

source: inundation/fit_interpolate/interpolate.py @ 2939

Last change on this file since 2939 was 2884, checked in by ole, 19 years ago
Replaced reference to f.T with f.get_time() and retested
File size: 26.2 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
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.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	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
238	n = point_coordinates.shape[0] #Nbr of data points
239
240	if verbose: print 'Number of datapoints: %d' %n
241	if verbose: print 'Number of basis functions: %d' %m
242
243	A = Sparse(n,m)
244
245	n = len(self.inside_poly_indices)
246	#Compute matrix elements for points inside the mesh
247	for i in self.inside_poly_indices:
248	#For each data_coordinate point
249	if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
250	x = point_coordinates[i]
251	element_found, sigma0, sigma1, sigma2, k = \
252	search_tree_of_vertices(self.root, self.mesh, x)
253	#Update interpolation matrix A if necessary
254	if element_found is True:
255	#Assign values to matrix A
256
257	j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
258	j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
259	j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
260
261	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
262	js = [j0,j1,j2]
263
264	for j in js:
265	A[i,j] = sigmas[j]
266	else:
267	msg = 'Could not find triangle for point', x
268	raise Exception(msg)
269	return A
270
271	class Interpolation_function:
272	"""Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
273	which is interpolated from time series defined at vertices of
274	triangular mesh (such as those stored in sww files)
275
276	Let m be the number of vertices, n the number of triangles
277	and p the number of timesteps.
278
279	Mandatory input
280	time: px1 array of monotonously increasing times (Float)
281	quantities: Dictionary of arrays or 1 array (Float)
282	The arrays must either have dimensions pxm or mx1.
283	The resulting function will be time dependent in
284	the former case while it will be constan with
285	respect to time in the latter case.
286
287	Optional input:
288	quantity_names: List of keys into the quantities dictionary
289	vertex_coordinates: mx2 array of coordinates (Float)
290	triangles: nx3 array of indices into vertex_coordinates (Int)
291	interpolation_points: Nx2 array of coordinates to be interpolated to
292	verbose: Level of reporting
293
294
295	The quantities returned by the callable object are specified by
296	the list quantities which must contain the names of the
297	quantities to be returned and also reflect the order, e.g. for
298	the shallow water wave equation, on would have
299	quantities = ['stage', 'xmomentum', 'ymomentum']
300
301	The parameter interpolation_points decides at which points interpolated
302	quantities are to be computed whenever object is called.
303	If None, return average value
304	"""
305
306
307
308	def __init__(self,
309	time,
310	quantities,
311	quantity_names = None,
312	vertex_coordinates = None,
313	triangles = None,
314	interpolation_points = None,
315	verbose = False):
316	"""Initialise object and build spatial interpolation if required
317	"""
318
319	from Numeric import array, zeros, Float, alltrue, concatenate,\
320	reshape, ArrayType
321
322
323	#from util import mean, ensure_numeric
324	from config import time_format
325	import types
326
327
328	#Check temporal info
329	time = ensure_numeric(time)
330	msg = 'Time must be a monotonuosly '
331	msg += 'increasing sequence %s' %time
332	assert alltrue(time[1:] - time[:-1] >= 0 ), msg
333
334
335	#Check if quantities is a single array only
336	if type(quantities) != types.DictType:
337	quantities = ensure_numeric(quantities)
338	quantity_names = ['Attribute']
339
340	#Make it a dictionary
341	quantities = {quantity_names[0]: quantities}
342
343
344	#Use keys if no names are specified
345	if quantity_names is None:
346	quantity_names = quantities.keys()
347
348
349	#Check spatial info
350	if vertex_coordinates is None:
351	self.spatial = False
352	else:
353	vertex_coordinates = ensure_numeric(vertex_coordinates)
354
355	assert triangles is not None, 'Triangles array must be specified'
356	triangles = ensure_numeric(triangles)
357	self.spatial = True
358
359
360	#Save for use with statistics
361	self.quantity_names = quantity_names
362	self.quantities = quantities
363	self.vertex_coordinates = vertex_coordinates
364	self.interpolation_points = interpolation_points
365	self.time = time[:] # Time assumed to be relative to starttime
366	self.index = 0 # Initial time index
367	self.precomputed_values = {}
368
369	#Precomputed spatial interpolation if requested
370	if interpolation_points is not None:
371	if self.spatial is False:
372	raise 'Triangles and vertex_coordinates must be specified'
373
374	try:
375	self.interpolation_points = ensure_numeric(interpolation_points)
376	except:
377	msg = 'Interpolation points must be an N x 2 Numeric array '+\
378	'or a list of points\n'
379	msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
380	'...')
381	raise msg
382
383
384	m = len(self.interpolation_points)
385	p = len(self.time)
386
387	for name in quantity_names:
388	self.precomputed_values[name] = zeros((p, m), Float)
389
390	#Build interpolator
391	interpol = Interpolate(vertex_coordinates,
392	triangles,
393	#point_coordinates = \
394	#self.interpolation_points,
395	#alpha = 0,
396	verbose = verbose)
397
398	if verbose: print 'Interpolate'
399	for i, t in enumerate(self.time):
400	#Interpolate quantities at this timestep
401	if verbose and i%((p+10)/10)==0:
402	print ' time step %d of %d' %(i, p)
403
404	for name in quantity_names:
405	if len(quantities[name].shape) == 2:
406	result = interpol.interpolate(quantities[name][i,:],
407	point_coordinates = \
408	self.interpolation_points)
409	else:
410	#Assume no time dependency
411	result = interpol.interpolate(quantities[name][:],
412	point_coordinates = \
413	self.interpolation_points)
414
415	self.precomputed_values[name][i, :] = result
416
417	#Report
418	if verbose:
419	print self.statistics()
420	#self.print_statistics()
421
422	else:
423	#Store quantitites as is
424	for name in quantity_names:
425	self.precomputed_values[name] = quantities[name]
426
427
428	def __repr__(self):
429	#return 'Interpolation function (spatio-temporal)'
430	return self.statistics()
431
432
433	def __call__(self, t, point_id = None, x = None, y = None):
434	"""Evaluate f(t), f(t, point_id) or f(t, x, y)
435
436	Inputs:
437	t: time - Model time. Must lie within existing timesteps
438	point_id: index of one of the preprocessed points.
439	x, y: Overrides location, point_id ignored
440
441	If spatial info is present and all of x,y,point_id
442	are None an exception is raised
443
444	If no spatial info is present, point_id and x,y arguments are ignored
445	making f a function of time only.
446
447
448	FIXME: point_id could also be a slice
449	FIXME: What if x and y are vectors?
450	FIXME: What about f(x,y) without t?
451	"""
452
453	from math import pi, cos, sin, sqrt
454	from Numeric import zeros, Float
455	from util import mean
456
457	if self.spatial is True:
458	if point_id is None:
459	if x is None or y is None:
460	msg = 'Either point_id or x and y must be specified'
461	raise Exception(msg)
462	else:
463	if self.interpolation_points is None:
464	msg = 'Interpolation_function must be instantiated ' +\
465	'with a list of interpolation points before parameter ' +\
466	'point_id can be used'
467	raise Exception(msg)
468
469	msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
470	msg += ' does not match model time: %.16f\n' %t
471	if t < self.time[0]: raise Exception(msg)
472	if t > self.time[-1]: raise Exception(msg)
473
474	oldindex = self.index #Time index
475
476	#Find current time slot
477	while t > self.time[self.index]: self.index += 1
478	while t < self.time[self.index]: self.index -= 1
479
480	if t == self.time[self.index]:
481	#Protect against case where t == T[-1] (last time)
482	# - also works in general when t == T[i]
483	ratio = 0
484	else:
485	#t is now between index and index+1
486	ratio = (t - self.time[self.index])/\
487	(self.time[self.index+1] - self.time[self.index])
488
489	#Compute interpolated values
490	q = zeros(len(self.quantity_names), Float)
491	#print "self.precomputed_values", self.precomputed_values
492	for i, name in enumerate(self.quantity_names):
493	Q = self.precomputed_values[name]
494
495	if self.spatial is False:
496	#If there is no spatial info
497	assert len(Q.shape) == 1
498
499	Q0 = Q[self.index]
500	if ratio > 0: Q1 = Q[self.index+1]
501
502	else:
503	if x is not None and y is not None:
504	#Interpolate to x, y
505
506	raise 'x,y interpolation not yet implemented'
507	else:
508	#Use precomputed point
509	Q0 = Q[self.index, point_id]
510	if ratio > 0:
511	Q1 = Q[self.index+1, point_id]
512
513	#Linear temporal interpolation
514	if ratio > 0:
515	if Q0 == INF and Q1 == INF:
516	q[i] = Q0
517	else:
518	q[i] = Q0 + ratio*(Q1 - Q0)
519	else:
520	q[i] = Q0
521
522
523	#Return vector of interpolated values
524	#if len(q) == 1:
525	# return q[0]
526	#else:
527	# return q
528
529
530	#Return vector of interpolated values
531	#FIXME:
532	if self.spatial is True:
533	return q
534	else:
535	#Replicate q according to x and y
536	#This is e.g used for Wind_stress
537	if x is None or y is None:
538	return q
539	else:
540	try:
541	N = len(x)
542	except:
543	return q
544	else:
545	from Numeric import ones, Float
546	#x is a vector - Create one constant column for each value
547	N = len(x)
548	assert len(y) == N, 'x and y must have same length'
549	res = []
550	for col in q:
551	res.append(col*ones(N, Float))
552
553	return res
554
555
556	def get_time(self):
557	"""Return model time as a vector of timesteps
558	"""
559	return self.time
560
561	def statistics(self):
562	"""Output statistics about interpolation_function
563	"""
564
565	vertex_coordinates = self.vertex_coordinates
566	interpolation_points = self.interpolation_points
567	quantity_names = self.quantity_names
568	quantities = self.quantities
569	precomputed_values = self.precomputed_values
570
571	x = vertex_coordinates[:,0]
572	y = vertex_coordinates[:,1]
573
574	str = '------------------------------------------------\n'
575	str += 'Interpolation_function (spatio-temporal) statistics:\n'
576	str += ' Extent:\n'
577	str += ' x in [%f, %f], len(x) == %d\n'\
578	%(min(x), max(x), len(x))
579	str += ' y in [%f, %f], len(y) == %d\n'\
580	%(min(y), max(y), len(y))
581	str += ' t in [%f, %f], len(t) == %d\n'\
582	%(min(self.time), max(self.time), len(self.time))
583	str += ' Quantities:\n'
584	for name in quantity_names:
585	q = quantities[name][:].flat
586	str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
587
588	if interpolation_points is not None:
589	str += ' Interpolation points (xi, eta):'\
590	' number of points == %d\n' %interpolation_points.shape[0]
591	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
592	max(interpolation_points[:,0]))
593	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
594	max(interpolation_points[:,1]))
595	str += ' Interpolated quantities (over all timesteps):\n'
596
597	for name in quantity_names:
598	q = precomputed_values[name][:].flat
599	str += ' %s at interpolation points in [%f, %f]\n'\
600	%(name, min(q), max(q))
601	str += '------------------------------------------------\n'
602
603	return str
604
605	def interpolate_sww(sww_file, time, interpolation_points,
606	quantity_names = None, verbose = False):
607	"""
608	obsolete.
609	use file_function in utils
610	"""
611	#open sww file
612	x, y, volumes, time, quantities = read_sww(sww_file)
613	print "x",x
614	print "y",y
615
616	print "time", time
617	print "quantities", quantities
618
619	#Add the x and y together
620	vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
621
622	#Will return the quantity values at the specified times and locations
623	interp = Interpolation_interface(
624	time,
625	quantities,
626	quantity_names = quantity_names,
627	vertex_coordinates = vertex_coordinates,
628	triangles = volumes,
629	interpolation_points = interpolation_points,
630	verbose = verbose)
631
632
633	def read_sww(file_name):
634	"""
635	obsolete - Nothing should be calling this
636
637	Read in an sww file.
638
639	Input;
640	file_name - the sww file
641
642	Output;
643	x - Vector of x values
644	y - Vector of y values
645	z - Vector of bed elevation
646	volumes - Array. Each row has 3 values, representing
647	the vertices that define the volume
648	time - Vector of the times where there is stage information
649	stage - array with respect to time and vertices (x,y)
650	"""
651
652	#FIXME Have this reader as part of data_manager?
653
654	from Scientific.IO.NetCDF import NetCDFFile
655	import tempfile
656	import sys
657	import os
658
659	#Check contents
660	#Get NetCDF
661
662	# see if the file is there. Throw a QUIET IO error if it isn't
663	# I don't think I could implement the above
664
665	#throws prints to screen if file not present
666	junk = tempfile.mktemp(".txt")
667	fd = open(junk,'w')
668	stdout = sys.stdout
669	sys.stdout = fd
670	fid = NetCDFFile(file_name, 'r')
671	sys.stdout = stdout
672	fd.close()
673	#clean up
674	os.remove(junk)
675
676	# Get the variables
677	x = fid.variables['x'][:]
678	y = fid.variables['y'][:]
679	volumes = fid.variables['volumes'][:]
680	time = fid.variables['time'][:]
681
682	keys = fid.variables.keys()
683	keys.remove("x")
684	keys.remove("y")
685	keys.remove("volumes")
686	keys.remove("time")
687	#Turn NetCDF objects into Numeric arrays
688	quantities = {}
689	for name in keys:
690	quantities[name] = fid.variables[name][:]
691
692	fid.close()
693	return x, y, volumes, time, quantities
694
695	#-------------------------------------------------------------
696	if __name__ == "__main__":
697	a = [0.0, 0.0]
698	b = [0.0, 2.0]
699	c = [2.0,0.0]
700	points = [a, b, c]
701	vertices = [ [1,0,2] ] #bac
702
703	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
704
705	interp = Interpolate(points, vertices) #, data)
706	A = interp._build_interpolation_matrix_A(data, verbose=True)
707	A = A.todense()
708	print "A",A
709	assert allclose(A, [[1./3, 1./3, 1./3]])
710	print "finished"

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

Download in other formats: