Context Navigation

source: inundation/fit_interpolate/interpolate.py @ 2691

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

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

Download in other formats: