Context Navigation

source: trunk/anuga_core/source/anuga/fit_interpolate/interpolate.py @ 8125

Last change on this file since 8125 was 8125, checked in by wilsonr, 13 years ago
Changes to address ticket 360.
File size: 48.4 KB

Line
1	"""Least squares interpolation.
2
3	These functions and classes calculate a value at a particular point on
4	the given mesh. It interpolates the values stored at the vertices of the
5	mesh.
6
7	For example, if you want to get the height of a terrain mesh at particular
8	point, you pass the point to an Interpolate class. The point will intersect
9	one of the triangles on the mesh, and the interpolated height will be an
10	intermediate value between the three vertices of that triangle.
11	This value is returned by the class.
12
13	Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
14	Geoscience Australia, 2004.
15
16	DESIGN ISSUES
17	* what variables should be global?
18	- if there are no global vars functions can be moved around alot easier
19
20	* The public interface to Interpolate
21	__init__
22	interpolate
23	interpolate_block
24
25	"""
26
27	import time
28	import os
29	import sys
30	from warnings import warn
31	from math import sqrt
32	from csv import writer, DictWriter
33
34	from anuga.caching.caching import cache
35	from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
36	from anuga.utilities.sparse import Sparse, Sparse_CSR
37	from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError
38	from anuga.coordinate_transforms.geo_reference import Geo_reference
39	from anuga.utilities.numerical_tools import ensure_numeric, NAN
40	from anuga.geospatial_data.geospatial_data import Geospatial_data
41	from anuga.geospatial_data.geospatial_data import ensure_absolute
42	from anuga.pmesh.mesh_quadtree import MeshQuadtree
43	from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
44	from anuga.abstract_2d_finite_volumes.file_function import file_function
45	from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
46	from anuga.geometry.polygon import interpolate_polyline, in_and_outside_polygon
47	import anuga.utilities.log as log
48
49	import numpy as num
50
51
52	# Interpolation specific exceptions
53
54	class Modeltime_too_late(Exception): pass
55	class Modeltime_too_early(Exception): pass
56
57
58	##
59	# @brief Interpolate vertex_values to interpolation points.
60	# @param vertex_coordinates List of coordinate pairs making a mesh.
61	# @param triangles Iterable of 3-tuples representing indices of mesh vertices.
62	# @param vertex_values Array of data at mesh vertices.
63	# @param interpolation_points Points to interpolate to.
64	# @param mesh_origin A geo_ref object or 3-tuples of UTMzone, easting, northing.
65	# @param start_blocking_len Block if # of points greater than this.
66	# @param use_cache If True, cache.
67	# @param verbose True if this function is to be verbose.
68	def interpolate(vertex_coordinates,
69	triangles,
70	vertex_values,
71	interpolation_points,
72	mesh_origin=None,
73	start_blocking_len=500000,
74	use_cache=False,
75	verbose=False,
76	output_centroids=False):
77	"""Interpolate vertex_values to interpolation points.
78
79	Inputs (mandatory):
80
81
82	vertex_coordinates: List of coordinate pairs [xi, eta] of
83	points constituting a mesh
84	(or an m x 2 numeric array or
85	a geospatial object)
86	Points may appear multiple times
87	(e.g. if vertices have discontinuities)
88
89	triangles: List of 3-tuples (or a numeric array) of
90	integers representing indices of all vertices
91	in the mesh.
92
93	vertex_values: Vector or array of data at the mesh vertices.
94	If array, interpolation will be done for each column as
95	per underlying matrix-matrix multiplication
96
97	interpolation_points: Interpolate mesh data to these positions.
98	List of coordinate pairs [x, y] of
99	data points or an nx2 numeric array or a
100	Geospatial_data object
101
102	Inputs (optional)
103
104	mesh_origin: A geo_reference object or 3-tuples consisting of
105	UTM zone, easting and northing.
106	If specified vertex coordinates are assumed to be
107	relative to their respective origins.
108
109	Note: Don't supply a vertex coords as a geospatial
110	object and a mesh origin, since geospatial has its
111	own mesh origin.
112
113	start_blocking_len: If the # of points is more or greater than this,
114	start blocking
115
116	use_cache: True or False
117
118
119	Output:
120
121	Interpolated values at specified point_coordinates
122
123	Note: This function is a simple shortcut for case where
124	interpolation matrix is unnecessary
125	Note: This function does not take blocking into account,
126	but allows caching.
127
128	"""
129
130	# FIXME(Ole): Probably obsolete since I is precomputed and
131	# interpolate_block caches
132
133	from anuga.caching import cache
134
135	# Create interpolation object with matrix
136	args = (ensure_numeric(vertex_coordinates, num.float),
137	ensure_numeric(triangles))
138	kwargs = {'mesh_origin': mesh_origin,
139	'verbose': verbose}
140
141	if use_cache is True:
142	if sys.platform != 'win32':
143	I = cache(Interpolate, args, kwargs, verbose=verbose)
144	else:
145	# Messy wrapping of Interpolate to deal with win32 error
146	def wrap_Interpolate(args,kwargs):
147	I = apply(Interpolate, args, kwargs)
148	return I
149	I = cache(wrap_Interpolate, (args, kwargs), {}, verbose=verbose)
150	else:
151	I = apply(Interpolate, args, kwargs)
152
153	# Call interpolate method with interpolation points
154	result = I.interpolate_block(vertex_values, interpolation_points,
155	use_cache=use_cache,
156	verbose=verbose,
157	output_centroids=output_centroids)
158
159	return result
160
161
162	##
163	# @brief
164	class Interpolate (FitInterpolate):
165
166	##
167	# @brief Build interpolation matrix.
168	# @param vertex_coordinates List of pairs [xi, eta] of points making a mesh.
169	# @param triangles List of 3-tuples of indices of all vertices in the mesh.
170	# @param mesh_origin A geo_ref object (UTM zone, easting and northing).
171	# @param verbose If True, this function is to be verbose.
172	# @param max_vertices_per_cell Split quadtree cell if vertices >= this.
173	def __init__(self,
174	vertex_coordinates,
175	triangles,
176	mesh_origin=None,
177	verbose=False):
178
179	""" Build interpolation matrix mapping from
180	function values at vertices to function values at data points
181
182	Inputs:
183	vertex_coordinates: List of coordinate pairs [xi, eta] of
184	points constituting a mesh (or an m x 2 numeric array or
185	a geospatial object)
186	Points may appear multiple times
187	(e.g. if vertices have discontinuities)
188
189	triangles: List of 3-tuples (or a numeric array) of
190	integers representing indices of all vertices in the mesh.
191
192	mesh_origin: A geo_reference object or 3-tuples consisting of
193	UTM zone, easting and northing.
194	If specified vertex coordinates are assumed to be
195	relative to their respective origins.
196
197	max_vertices_per_cell: Number of vertices in a quad tree cell
198	at which the cell is split into 4.
199
200	Note: Don't supply a vertex coords as a geospatial object and
201	a mesh origin, since geospatial has its own mesh origin.
202	"""
203
204	# FIXME (Ole): Need an input check
205
206	FitInterpolate.__init__(self,
207	vertex_coordinates=vertex_coordinates,
208	triangles=triangles,
209	mesh_origin=mesh_origin,
210	verbose=verbose)
211
212	# Initialise variables
213	self._A_can_be_reused = False # FIXME (Ole): Probably obsolete
214	self._point_coordinates = None # FIXME (Ole): Probably obsolete
215	self.interpolation_matrices = {} # Store precomputed matrices
216
217
218	##
219	# @brief Interpolate mesh data f to determine values, z, at points.
220	# @param f Data on the mesh vertices.
221	# @param point_coordinates Interpolate mesh data to these positions.
222	# @param start_blocking_len Block if # points >= this.
223	# @param verbose True if this function is to be verbose.
224	# FIXME: What is a good start_blocking_len value?
225	def interpolate(self,
226	f,
227	point_coordinates=None,
228	start_blocking_len=500000,
229	verbose=False,
230	output_centroids=False):
231	"""Interpolate mesh data f to determine values, z, at points.
232
233	f is the data on the mesh vertices.
234
235	The mesh values representing a smooth surface are
236	assumed to be specified in f.
237
238	Inputs:
239	f: Vector or array of data at the mesh vertices.
240	If f is an array, interpolation will be done for each column as
241	per underlying matrix-matrix multiplication
242
243	point_coordinates: Interpolate mesh data to these positions.
244	List of coordinate pairs [x, y] of
245	data points or an nx2 numeric array or a Geospatial_data object
246
247	If point_coordinates is absent, the points inputted last time
248	this method was called are used, if possible.
249
250	start_blocking_len: If the # of points is more or greater than this,
251	start blocking
252
253	Output:
254	Interpolated values at inputted points (z).
255	"""
256
257	# FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
258	# method is called even if interpolation points are unchanged.
259	# This really should use some kind of caching in cases where
260	# interpolation points are reused.
261	#
262	# This has now been addressed through an attempt in interpolate_block
263
264	if verbose: log.critical('Build intepolation object')
265	if isinstance(point_coordinates, Geospatial_data):
266	point_coordinates = point_coordinates.get_data_points(absolute=True)
267
268	# Can I interpolate, based on previous point_coordinates?
269	if point_coordinates is None:
270	if self._A_can_be_reused is True \
271	and len(self._point_coordinates) < start_blocking_len:
272	z = self._get_point_data_z(f, verbose=verbose)
273	elif self._point_coordinates is not None:
274	# if verbose, give warning
275	if verbose:
276	log.critical('WARNING: Recalculating A matrix, '
277	'due to blocking.')
278	point_coordinates = self._point_coordinates
279	else:
280	# There are no good point_coordinates. import sys; sys.exit()
281	msg = 'ERROR (interpolate.py): No point_coordinates inputted'
282	raise Exception(msg)
283
284	if point_coordinates is not None:
285	self._point_coordinates = point_coordinates
286	if len(point_coordinates) < start_blocking_len \
287	or start_blocking_len == 0:
288	self._A_can_be_reused = True
289	z = self.interpolate_block(f, point_coordinates,
290	verbose=verbose, output_centroids=output_centroids)
291	else:
292	# Handle blocking
293	self._A_can_be_reused = False
294	start = 0
295	# creating a dummy array to concatenate to.
296
297	f = ensure_numeric(f, num.float)
298	if len(f.shape) > 1:
299	z = num.zeros((0, f.shape[1]), num.int) #array default#
300	else:
301	z = num.zeros((0,), num.int) #array default#
302
303	for end in range(start_blocking_len,
304	len(point_coordinates),
305	start_blocking_len):
306	t = self.interpolate_block(f, point_coordinates[start:end],
307	verbose=verbose, output_centroids=output_centroids)
308	z = num.concatenate((z, t), axis=0) #??default#
309	start = end
310
311	end = len(point_coordinates)
312	t = self.interpolate_block(f, point_coordinates[start:end],
313	verbose=verbose, output_centroids=output_centroids)
314	z = num.concatenate((z, t), axis=0) #??default#
315	return z
316
317
318	##
319	# @brief Interpolate a block of vertices
320	# @param f Array of arbitrary data to be interpolated
321	# @param point_coordinates List of vertices to intersect with the mesh
322	# @param use_cache True if caching should be used to accelerate the calculations
323	# @param verbose True if this function is verbose.
324	# @return interpolated f
325	def interpolate_block(self, f, point_coordinates,
326	use_cache=False, verbose=False, output_centroids=False):
327	"""
328	Call this if you want to control the blocking or make sure blocking
329	doesn't occur.
330
331	Return the point data, z.
332
333	See interpolate for doc info.
334	"""
335
336	# FIXME (Ole): I reckon we should change the interface so that
337	# the user can specify the interpolation matrix instead of the
338	# interpolation points to save time.
339
340	if isinstance(point_coordinates, Geospatial_data):
341	point_coordinates = point_coordinates.get_data_points(absolute=True)
342
343	# Convert lists to numeric arrays if necessary
344	point_coordinates = ensure_numeric(point_coordinates, num.float)
345	f = ensure_numeric(f, num.float)
346
347	from anuga.caching import myhash
348	import sys
349
350	if use_cache is True:
351	if sys.platform != 'win32':
352	# FIXME (Ole): (Why doesn't this work on windoze?)
353	# Still absolutely fails on Win 24 Oct 2008
354
355	X = cache(self._build_interpolation_matrix_A,
356	args=(point_coordinates, output_centroids),
357	kwargs={'verbose': verbose},
358	verbose=verbose)
359	else:
360	# FIXME
361	# Hash point_coordinates to memory location, reuse if possible
362	# This will work on Linux as well if we want to use it there.
363	key = myhash(point_coordinates)
364
365	reuse_A = False
366
367	if self.interpolation_matrices.has_key(key):
368	X, stored_points = self.interpolation_matrices[key]
369	if num.alltrue(stored_points == point_coordinates):
370	reuse_A = True # Reuse interpolation matrix
371
372	if reuse_A is False:
373	X = self._build_interpolation_matrix_A(point_coordinates,
374	output_centroids,
375	verbose=verbose)
376	self.interpolation_matrices[key] = (X, point_coordinates)
377	else:
378	X = self._build_interpolation_matrix_A(point_coordinates, output_centroids,
379	verbose=verbose)
380
381	# Unpack result
382	self._A, self.inside_poly_indices, self.outside_poly_indices, self.centroids = X
383
384	# Check that input dimensions are compatible
385	msg = 'Two columns must be specified in point coordinates. ' \
386	'I got shape=%s' % (str(point_coordinates.shape))
387	assert point_coordinates.shape[1] == 2, msg
388
389	msg = 'The number of rows in matrix A must be the same as the '
390	msg += 'number of points supplied.'
391	msg += ' I got %d points and %d matrix rows.' \
392	% (point_coordinates.shape[0], self._A.shape[0])
393	assert point_coordinates.shape[0] == self._A.shape[0], msg
394
395	msg = 'The number of columns in matrix A must be the same as the '
396	msg += 'number of mesh vertices.'
397	msg += ' I got %d vertices and %d matrix columns.' \
398	% (f.shape[0], self._A.shape[1])
399	assert self._A.shape[1] == f.shape[0], msg
400
401	# Compute Matrix vector product and return
402	return self._get_point_data_z(f)
403
404
405	##
406	# @brief Get interpolated data at given points.
407	# Applies a transform to all points to calculate the
408	# interpolated values. Points outside the mesh are returned as NaN.
409	# @note self._A matrix must be valid
410	# @param f Array of arbitrary data
411	# @param verbose True if this function is to be verbose.
412	# @return f transformed by interpolation matrix (f')
413	def _get_point_data_z(self, f, verbose=False):
414	"""
415	Return the point data, z.
416
417	Precondition: The _A matrix has been created
418	"""
419
420	z = self._A * f
421
422	# Taking into account points outside the mesh.
423	for i in self.outside_poly_indices:
424	z[i] = NAN
425	return z
426
427
428	##
429	# @brief Build NxM interpolation matrix.
430	# @param point_coordinates Points to sample at
431	# @param output_centroids set to True to always sample from the centre
432	# of the intersected triangle, instead of the intersection
433	# point.
434	# @param verbose True if this function is to be verbose.
435	# @return Interpolation matrix A, plus lists of the points inside and outside the mesh
436	# and the list of centroids, if requested.
437	def _build_interpolation_matrix_A(self,
438	point_coordinates,
439	output_centroids=False,
440	verbose=False):
441	"""Build n x m interpolation matrix, where
442	n is the number of data points and
443	m is the number of basis functions phi_k (one per vertex)
444
445	This algorithm uses a quad tree data structure for fast binning
446	of data points
447	origin is a 3-tuple consisting of UTM zone, easting and northing.
448	If specified coordinates are assumed to be relative to this origin.
449
450	This one will override any data_origin that may be specified in
451	instance interpolation
452
453	Preconditions:
454	Point_coordindates and mesh vertices have the same origin.
455	"""
456
457	if verbose: log.critical('Building interpolation matrix')
458
459	# Convert point_coordinates to numeric arrays, in case it was a list.
460	point_coordinates = ensure_numeric(point_coordinates, num.float)
461
462	if verbose: log.critical('Getting indices inside mesh boundary')
463
464	# Quick test against boundary, but will not deal with holes in the mesh
465	inside_boundary_indices, outside_poly_indices = \
466	in_and_outside_polygon(point_coordinates,
467	self.mesh.get_boundary_polygon(),
468	closed=True, verbose=verbose)
469
470	# Build n x m interpolation matrix
471	if verbose and len(outside_poly_indices) > 0:
472	log.critical('WARNING: Points outside mesh boundary.')
473
474	# Since you can block, throw a warning, not an error.
475	if verbose and 0 == len(inside_boundary_indices):
476	log.critical('WARNING: No points within the mesh!')
477
478	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
479	n = point_coordinates.shape[0] # Nbr of data points
480
481	if verbose: log.critical('Number of datapoints: %d' % n)
482	if verbose: log.critical('Number of basis functions: %d' % m)
483
484	A = Sparse(n,m)
485
486	n = len(inside_boundary_indices)
487
488	centroids = []
489	inside_poly_indices = []
490
491	# Compute matrix elements for points inside the mesh
492	if verbose: log.critical('Building interpolation matrix from %d points'
493	% n)
494
495	for d, i in enumerate(inside_boundary_indices):
496	# For each data_coordinate point
497	if verbose and d%((n+10)/10)==0: log.critical('Doing %d of %d'
498	%(d, n))
499
500	x = point_coordinates[i]
501	element_found, sigma0, sigma1, sigma2, k = self.root.search_fast(x)
502
503	# Update interpolation matrix A if necessary
504	if element_found is True:
505
506	if verbose:
507	print 'Point is within mesh:', d, i
508
509	inside_poly_indices.append(i)
510
511	# Assign values to matrix A
512	j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0
513	j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
514	j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
515	js = [j0, j1, j2]
516
517	if output_centroids is False:
518	# Weight each vertex according to its distance from x
519	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
520	for j in js:
521	A[i, j] = sigmas[j]
522	else:
523	# If centroids are needed, weight all 3 vertices equally
524	for j in js:
525	A[i, j] = 1.0/3.0
526	centroids.append(self.mesh.centroid_coordinates[k])
527	else:
528	if verbose:
529	log.critical('Mesh has a hole - moving this point to outside list')
530
531	# This is a numpy arrays, so we need to do a slow transfer
532	outside_poly_indices = num.append(outside_poly_indices, [i], axis=0)
533
534	return A, inside_poly_indices, outside_poly_indices, centroids
535
536
537
538
539
540
541	##
542	# @brief ??
543	# @param vertices ??
544	# @param vertex_attributes ??
545	# @param triangles ??
546	# @param points ??
547	# @param max_points_per_cell ??
548	# @param start_blocking_len ??
549	# @param mesh_origin ??
550	def benchmark_interpolate(vertices,
551	vertex_attributes,
552	triangles, points,
553	max_points_per_cell=None,
554	start_blocking_len=500000,
555	mesh_origin=None):
556	"""
557	points: Interpolate mesh data to these positions.
558	List of coordinate pairs [x, y] of
559	data points or an nx2 numeric array or a Geospatial_data object
560
561	No test for this yet.
562	Note, this has no time the input data has no time dimension. Which is
563	different from most of the data we interpolate, eg sww info.
564
565	Output:
566	Interpolated values at inputted points.
567	"""
568
569	interp = Interpolate(vertices,
570	triangles,
571	max_vertices_per_cell=max_points_per_cell,
572	mesh_origin=mesh_origin)
573
574	calc = interp.interpolate(vertex_attributes,
575	points,
576	start_blocking_len=start_blocking_len)
577
578
579	##
580	# @brief Interpolate quantities at given locations (from .SWW file).
581	# @param sww_file Input .SWW file.
582	# @param points A list of the 'gauges' x,y location.
583	# @param depth_file The name of the output depth file.
584	# @param velocity_x_file Name of the output x velocity file.
585	# @param velocity_y_file Name of the output y velocity file.
586	# @param stage_file Name of the output stage file.
587	# @param froude_file
588	# @param time_thinning Time thinning step to use.
589	# @param verbose True if this function is to be verbose.
590	# @param use_cache True if we are caching.
591	def interpolate_sww2csv(sww_file,
592	points,
593	depth_file,
594	velocity_x_file,
595	velocity_y_file,
596	stage_file=None,
597	froude_file=None,
598	time_thinning=1,
599	verbose=True,
600	use_cache = True):
601	"""
602	Interpolate the quantities at a given set of locations, given
603	an sww file.
604	The results are written to csv files.
605
606	sww_file is the input sww file.
607	points is a list of the 'gauges' x,y location.
608	depth_file is the name of the output depth file
609	velocity_x_file is the name of the output x velocity file.
610	velocity_y_file is the name of the output y velocity file.
611	stage_file is the name of the output stage file.
612
613	In the csv files columns represents the gauges and each row is a
614	time slice.
615
616	Time_thinning_number controls how many timesteps to use. Only
617	timesteps with index%time_thinning_number == 0 will used, or
618	in other words a value of 3, say, will cause the algorithm to
619	use every third time step.
620
621	In the future let points be a points file.
622	And let the user choose the quantities.
623
624	This is currently quite specific.
625	If it is need to be more general, change things.
626	"""
627
628	quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
629	points = ensure_absolute(points)
630	point_count = len(points)
631	callable_sww = file_function(sww_file,
632	quantities=quantities,
633	interpolation_points=points,
634	verbose=verbose,
635	time_thinning=time_thinning,
636	use_cache=use_cache)
637
638	depth_writer = writer(file(depth_file, "wb"))
639	velocity_x_writer = writer(file(velocity_x_file, "wb"))
640	velocity_y_writer = writer(file(velocity_y_file, "wb"))
641	if stage_file is not None:
642	stage_writer = writer(file(stage_file, "wb"))
643	if froude_file is not None:
644	froude_writer = writer(file(froude_file, "wb"))
645
646	# Write heading
647	heading = [str(x[0])+ ':' + str(x[1]) for x in points]
648	heading.insert(0, "time")
649	depth_writer.writerow(heading)
650	velocity_x_writer.writerow(heading)
651	velocity_y_writer.writerow(heading)
652	if stage_file is not None:
653	stage_writer.writerow(heading)
654	if froude_file is not None:
655	froude_writer.writerow(heading)
656
657	for time in callable_sww.get_time():
658	depths = [time]
659	velocity_xs = [time]
660	velocity_ys = [time]
661	if stage_file is not None:
662	stages = [time]
663	if froude_file is not None:
664	froudes = [time]
665	for point_i, point in enumerate(points):
666	quantities = callable_sww(time,point_i)
667
668	w = quantities[0]
669	z = quantities[1]
670	momentum_x = quantities[2]
671	momentum_y = quantities[3]
672	depth = w - z
673
674	if w == NAN or z == NAN or momentum_x == NAN:
675	velocity_x = NAN
676	else:
677	if depth > 1.e-30: # use epsilon
678	velocity_x = momentum_x / depth #Absolute velocity
679	else:
680	velocity_x = 0
681
682	if w == NAN or z == NAN or momentum_y == NAN:
683	velocity_y = NAN
684	else:
685	if depth > 1.e-30: # use epsilon
686	velocity_y = momentum_y / depth #Absolute velocity
687	else:
688	velocity_y = 0
689
690	if depth < 1.e-30: # use epsilon
691	froude = NAN
692	else:
693	froude = sqrt(velocity_xvelocity_x + velocity_yvelocity_y)/ \
694	sqrt(depth * 9.8066) # gravity m/s/s
695
696	depths.append(depth)
697	velocity_xs.append(velocity_x)
698	velocity_ys.append(velocity_y)
699
700	if stage_file is not None:
701	stages.append(w)
702	if froude_file is not None:
703	froudes.append(froude)
704
705	depth_writer.writerow(depths)
706	velocity_x_writer.writerow(velocity_xs)
707	velocity_y_writer.writerow(velocity_ys)
708
709	if stage_file is not None:
710	stage_writer.writerow(stages)
711	if froude_file is not None:
712	froude_writer.writerow(froudes)
713
714
715	##
716	# @brief
717	class Interpolation_function:
718	"""Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
719	which is interpolated from time series defined at vertices of
720	triangular mesh (such as those stored in sww files)
721
722	Let m be the number of vertices, n the number of triangles
723	and p the number of timesteps.
724	Also, let N be the number of interpolation points.
725
726	Mandatory input
727	time: px1 array of monotonously increasing times (float)
728	quantities: Dictionary of arrays or 1 array (float)
729	The arrays must either have dimensions pxm or mx1.
730	The resulting function will be time dependent in
731	the former case while it will be constant with
732	respect to time in the latter case.
733
734	Optional input:
735	quantity_names: List of keys into the quantities dictionary for
736	imposing a particular order on the output vector.
737	vertex_coordinates: mx2 array of coordinates (float)
738	triangles: nx3 array of indices into vertex_coordinates (int)
739	interpolation_points: Nx2 array of coordinates to be interpolated to
740	verbose: Level of reporting
741
742	The quantities returned by the callable object are specified by
743	the list quantities which must contain the names of the
744	quantities to be returned and also reflect the order, e.g. for
745	the shallow water wave equation, on would have
746	quantities = ['stage', 'xmomentum', 'ymomentum']
747
748	The parameter interpolation_points decides at which points interpolated
749	quantities are to be computed whenever object is called.
750	If None, return average value
751
752	FIXME (Ole): Need to allow vertex coordinates and interpolation points to
753	be geospatial data objects
754
755	(FIXME (Ole): This comment should be removed)
756	Time assumed to be relative to starttime
757	All coordinates assume origin of (0,0) - e.g. georeferencing must be
758	taken care of outside this function
759	"""
760
761	##
762	# @brief ??
763	# @param time ??
764	# @param quantities ??
765	# @param quantity_names ??
766	# @param vertex_coordinates ??
767	# @param triangles ??
768	# @param interpolation_points ??
769	# @param time_thinning ??
770	# @param verbose ??
771	# @param gauge_neighbour_id ??
772	def __init__(self,
773	time,
774	quantities,
775	quantity_names=None,
776	vertex_coordinates=None,
777	triangles=None,
778	interpolation_points=None,
779	time_thinning=1,
780	verbose=False,
781	gauge_neighbour_id=None,
782	output_centroids=False):
783	"""Initialise object and build spatial interpolation if required
784
785	Time_thinning_number controls how many timesteps to use. Only timesteps
786	with index%time_thinning_number == 0 will used, or in other words a
787	value of 3, say, will cause the algorithm to use every third time step.
788	"""
789
790	from anuga.config import time_format
791
792	if verbose is True:
793	log.critical('Interpolation_function: input checks')
794
795	# Check temporal info
796	time = ensure_numeric(time)
797
798	if not num.alltrue(time[1:] - time[:-1] >= 0):
799	# This message is time consuming to form due to the conversion of
800	msg = 'Time must be a monotonuosly increasing sequence %s' % time
801	raise Exception(msg)
802
803	# Check if quantities is a single array only
804	if not isinstance(quantities, dict):
805	quantities = ensure_numeric(quantities)
806	quantity_names = ['Attribute']
807
808	# Make it a dictionary
809	quantities = {quantity_names[0]: quantities}
810
811	# Use keys if no names are specified
812	if quantity_names is None:
813	quantity_names = quantities.keys()
814
815	# Check spatial info
816	if vertex_coordinates is None:
817	self.spatial = False
818	else:
819	# FIXME (Ole): Try ensure_numeric here -
820	# this function knows nothing about georefering.
821	vertex_coordinates = ensure_absolute(vertex_coordinates)
822
823	if triangles is not None:
824	triangles = ensure_numeric(triangles)
825	self.spatial = True
826
827	if verbose is True:
828	log.critical('Interpolation_function: thinning by %d'
829	% time_thinning)
830
831
832	# Thin timesteps if needed
833	# Note array() is used to make the thinned arrays contiguous in memory
834	self.time = num.array(time[::time_thinning])
835	for name in quantity_names:
836	if len(quantities[name].shape) == 2:
837	quantities[name] = num.array(quantities[name][::time_thinning,:])
838
839	if verbose is True:
840	log.critical('Interpolation_function: precomputing')
841
842	# Save for use with statistics
843	self.quantities_range = {}
844	for name in quantity_names:
845	q = quantities[name][:].flatten()
846	self.quantities_range[name] = [min(q), max(q)]
847
848	self.quantity_names = quantity_names
849	self.vertex_coordinates = vertex_coordinates
850	self.interpolation_points = interpolation_points
851
852	self.index = 0 # Initial time index
853	self.precomputed_values = {}
854	self.centroids = []
855
856	# Precomputed spatial interpolation if requested
857	if interpolation_points is not None:
858	#no longer true. sts files have spatial = True but
859	#if self.spatial is False:
860	# raise Exception('Triangles and vertex_coordinates must be specified')
861	#
862	try:
863	self.interpolation_points = \
864	interpolation_points = ensure_numeric(interpolation_points)
865	except:
866	msg = 'Interpolation points must be an N x 2 numeric array ' \
867	'or a list of points\n'
868	msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...')
869	raise Exception(msg)
870
871	# Ensure 'mesh_boundary_polygon' is defined
872	mesh_boundary_polygon = None
873
874	if triangles is not None and vertex_coordinates is not None:
875	# Check that all interpolation points fall within
876	# mesh boundary as defined by triangles and vertex_coordinates.
877	from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
878	from anuga.geometry.polygon import outside_polygon
879
880	# Create temporary mesh object from mesh info passed
881	# into this function.
882	mesh = Mesh(vertex_coordinates, triangles)
883	mesh_boundary_polygon = mesh.get_boundary_polygon()
884
885	indices = outside_polygon(interpolation_points,
886	mesh_boundary_polygon)
887
888	# Record result
889	#self.mesh_boundary_polygon = mesh_boundary_polygon
890	self.indices_outside_mesh = indices
891
892	# Report
893	if len(indices) > 0:
894	msg = 'Interpolation points in Interpolation function fall '
895	msg += 'outside specified mesh. Offending points:\n'
896	out_interp_pts = []
897	for i in indices:
898	msg += '%d: %s\n' % (i, interpolation_points[i])
899	out_interp_pts.append(
900	ensure_numeric(interpolation_points[i]))
901
902	if verbose is True:
903	import sys
904	from anuga.geometry.polygon import plot_polygons
905	title = ('Interpolation points fall '
906	'outside specified mesh')
907	plot_polygons([mesh_boundary_polygon,
908	interpolation_points,
909	out_interp_pts],
910	['line', 'point', 'outside'],
911	figname='points_boundary_out',
912	label=title)
913
914	# Joaquim Luis suggested this as an Exception, so
915	# that the user can now what the problem is rather than
916	# looking for NaN's. However, NANs are handy as they can
917	# be ignored leaving good points for continued processing.
918	if verbose:
919	log.critical(msg)
920	#raise Exception(msg)
921
922	elif triangles is None and vertex_coordinates is not None: #jj
923	#Dealing with sts file
924	pass
925	else:
926	raise Exception('Sww file function requires both triangles and '
927	'vertex_coordinates. sts file file function '
928	'requires the latter.')
929
930	# Plot boundary and interpolation points,
931	# but only if if 'mesh_boundary_polygon' has data.
932	if verbose is True and mesh_boundary_polygon is not None:
933	import sys
934	if sys.platform == 'win32':
935	from anuga.geometry.polygon import plot_polygons
936	title = ('Interpolation function: '
937	'Polygon and interpolation points')
938	plot_polygons([mesh_boundary_polygon,
939	interpolation_points],
940	['line', 'point'],
941	figname='points_boundary',
942	label=title)
943
944	m = len(self.interpolation_points)
945	p = len(self.time)
946
947	for name in quantity_names:
948	self.precomputed_values[name] = num.zeros((p, m), num.float)
949
950	if verbose is True:
951	log.critical('Build interpolator')
952
953
954	# Build interpolator
955	if triangles is not None and vertex_coordinates is not None:
956	if verbose:
957	msg = 'Building interpolation matrix from source mesh '
958	msg += '(%d vertices, %d triangles)' \
959	% (vertex_coordinates.shape[0],
960	triangles.shape[0])
961	log.critical(msg)
962
963	# This one is no longer needed for STS files
964	interpol = Interpolate(vertex_coordinates,
965	triangles,
966	verbose=verbose)
967
968	elif triangles is None and vertex_coordinates is not None:
969	if verbose:
970	log.critical('Interpolation from STS file')
971
972
973
974	if verbose:
975	log.critical('Interpolating (%d interpolation points, %d timesteps).'
976	% (self.interpolation_points.shape[0], self.time.shape[0]))
977
978	if time_thinning > 1:
979	log.critical('Timesteps were thinned by a factor of %d'
980	% time_thinning)
981	else:
982	log.critical()
983
984	for i, t in enumerate(self.time):
985	# Interpolate quantities at this timestep
986	if verbose and i%((p+10)/10) == 0:
987	log.critical(' time step %d of %d' % (i, p))
988
989	for name in quantity_names:
990	if len(quantities[name].shape) == 2:
991	Q = quantities[name][i,:] # Quantities at timestep i
992	else:
993	Q = quantities[name][:] # No time dependency
994
995	if verbose and i%((p+10)/10) == 0:
996	log.critical(' quantity %s, size=%d' % (name, len(Q)))
997
998	# Interpolate
999	if triangles is not None and vertex_coordinates is not None:
1000	result = interpol.interpolate(Q,
1001	point_coordinates=\
1002	self.interpolation_points,
1003	verbose=False,
1004	output_centroids=output_centroids)
1005	self.centroids = interpol.centroids
1006	elif triangles is None and vertex_coordinates is not None:
1007	result = interpolate_polyline(Q,
1008	vertex_coordinates,
1009	gauge_neighbour_id,
1010	interpolation_points=\
1011	self.interpolation_points)
1012
1013	#assert len(result), len(interpolation_points)
1014	self.precomputed_values[name][i, :] = result
1015
1016	# Report
1017	if verbose:
1018	log.critical(self.statistics())
1019	else:
1020	# Store quantitites as is
1021	for name in quantity_names:
1022	self.precomputed_values[name] = quantities[name]
1023
1024	##
1025	# @brief Override object representation method.
1026	def __repr__(self):
1027	# return 'Interpolation function (spatio-temporal)'
1028	return self.statistics()
1029
1030	##
1031	# @brief Evaluate interpolation function
1032	# @param t Model time - must lie within existing timesteps.
1033	# @param point_id Index of one of the preprocessed points.
1034	# @param x ??
1035	# @param y ??
1036	# @return ??
1037	def __call__(self, t, point_id=None, x=None, y=None):
1038	"""Evaluate f(t) or f(t, point_id)
1039
1040	Inputs:
1041	t: time - Model time. Must lie within existing timesteps
1042	point_id: index of one of the preprocessed points.
1043
1044	If spatial info is present and all of point_id
1045	are None an exception is raised
1046
1047	If no spatial info is present, point_id arguments are ignored
1048	making f a function of time only.
1049
1050	FIXME: f(t, x, y) x, y could overrided location, point_id ignored
1051	FIXME: point_id could also be a slice
1052	FIXME: What if x and y are vectors?
1053	FIXME: What about f(x,y) without t?
1054	"""
1055
1056	from math import pi, cos, sin, sqrt
1057
1058	if self.spatial is True:
1059	if point_id is None:
1060	if x is None or y is None:
1061	msg = 'Either point_id or x and y must be specified'
1062	raise Exception(msg)
1063	else:
1064	if self.interpolation_points is None:
1065	msg = 'Interpolation_function must be instantiated ' + \
1066	'with a list of interpolation points before ' + \
1067	'parameter point_id can be used'
1068	raise Exception(msg)
1069
1070	msg = 'Time interval [%.16f:%.16f]' % (self.time[0], self.time[-1])
1071	msg += ' does not match model time: %.16f\n' % t
1072	if t < self.time[0]: raise Modeltime_too_early(msg)
1073	if t > self.time[-1]: raise Modeltime_too_late(msg)
1074
1075	oldindex = self.index #Time index
1076
1077	# Find current time slot
1078	while t > self.time[self.index]: self.index += 1
1079	while t < self.time[self.index]: self.index -= 1
1080
1081	if t == self.time[self.index]:
1082	# Protect against case where t == T[-1] (last time)
1083	# - also works in general when t == T[i]
1084	ratio = 0
1085	else:
1086	# t is now between index and index+1
1087	ratio = ((t - self.time[self.index]) /
1088	(self.time[self.index+1] - self.time[self.index]))
1089
1090	# Compute interpolated values
1091	q = num.zeros(len(self.quantity_names), num.float)
1092	for i, name in enumerate(self.quantity_names):
1093	Q = self.precomputed_values[name]
1094
1095	if self.spatial is False:
1096	# If there is no spatial info
1097	assert len(Q.shape) == 1
1098
1099	Q0 = Q[self.index]
1100	if ratio > 0: Q1 = Q[self.index+1]
1101	else:
1102	if x is not None and y is not None:
1103	# Interpolate to x, y
1104	raise Exception('x,y interpolation not yet implemented')
1105	else:
1106	# Use precomputed point
1107	Q0 = Q[self.index, point_id]
1108	if ratio > 0:
1109	Q1 = Q[self.index+1, point_id]
1110
1111	# Linear temporal interpolation
1112	if ratio > 0:
1113	if Q0 == NAN and Q1 == NAN:
1114	q[i] = Q0
1115	else:
1116	q[i] = Q0 + ratio*(Q1 - Q0)
1117	else:
1118	q[i] = Q0
1119
1120	# Return vector of interpolated values
1121	# FIXME:
1122	if self.spatial is True:
1123	return q
1124	else:
1125	# Replicate q according to x and y
1126	# This is e.g used for Wind_stress
1127	if x is None or y is None:
1128	return q
1129	else:
1130	try:
1131	N = len(x)
1132	except:
1133	return q
1134	else:
1135	# x is a vector - Create one constant column for each value
1136	N = len(x)
1137	assert len(y) == N, 'x and y must have same length'
1138	res = []
1139	for col in q:
1140	res.append(col*num.ones(N, num.float))
1141
1142	return res
1143
1144	##
1145	# @brief Return model time as a vector of timesteps.
1146	def get_time(self):
1147	"""Return model time as a vector of timesteps
1148	"""
1149	return self.time
1150
1151	##
1152	# @brief Output statistics about interpolation_function.
1153	# @return The statistics string.
1154	def statistics(self):
1155	"""Output statistics about interpolation_function
1156	"""
1157
1158	vertex_coordinates = self.vertex_coordinates
1159	interpolation_points = self.interpolation_points
1160	quantity_names = self.quantity_names
1161	#quantities = self.quantities
1162	precomputed_values = self.precomputed_values
1163
1164	x = vertex_coordinates[:,0]
1165	y = vertex_coordinates[:,1]
1166
1167	str = '------------------------------------------------\n'
1168	str += 'Interpolation_function (spatio-temporal) statistics:\n'
1169	str += ' Extent:\n'
1170	str += ' x in [%f, %f], len(x) == %d\n'\
1171	%(min(x), max(x), len(x))
1172	str += ' y in [%f, %f], len(y) == %d\n'\
1173	%(min(y), max(y), len(y))
1174	str += ' t in [%f, %f], len(t) == %d\n'\
1175	%(min(self.time), max(self.time), len(self.time))
1176	str += ' Quantities:\n'
1177	for name in quantity_names:
1178	minq, maxq = self.quantities_range[name]
1179	str += ' %s in [%f, %f]\n' %(name, minq, maxq)
1180	#q = quantities[name][:].flatten()
1181	#str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
1182
1183	if interpolation_points is not None:
1184	str += ' Interpolation points (xi, eta):'\
1185	' number of points == %d\n' %interpolation_points.shape[0]
1186	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
1187	max(interpolation_points[:,0]))
1188	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
1189	max(interpolation_points[:,1]))
1190	str += ' Interpolated quantities (over all timesteps):\n'
1191
1192	for name in quantity_names:
1193	q = precomputed_values[name][:].flatten()
1194	str += ' %s at interpolation points in [%f, %f]\n'\
1195	%(name, min(q), max(q))
1196	str += '------------------------------------------------\n'
1197
1198	return str
1199
1200
1201	##
1202	# @brief ??
1203	# @param sww_file ??
1204	# @param time ??
1205	# @param interpolation_points ??
1206	# @param quantity_names ??
1207	# @param verbose ??
1208	# @note Obsolete. Use file_function() in utils.
1209	def interpolate_sww(sww_file, time, interpolation_points,
1210	quantity_names=None, verbose=False):
1211	"""
1212	obsolete.
1213	use file_function in utils
1214	"""
1215
1216	#open sww file
1217	x, y, volumes, time, quantities = read_sww(sww_file)
1218	log.critical("x=%s" % str(x))
1219	log.critical("y=%s" % str(y))
1220
1221	log.critical("time=%s" % str(time))
1222	log.critical("quantities=%s" % str(quantities))
1223
1224	#Add the x and y together
1225	vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),
1226	axis=1)
1227
1228	#Will return the quantity values at the specified times and locations
1229	interp = Interpolation_interface(time,
1230	quantities,
1231	quantity_names=quantity_names,
1232	vertex_coordinates=vertex_coordinates,
1233	triangles=volumes,
1234	interpolation_points=interpolation_points,
1235	verbose=verbose)
1236
1237

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

Download in other formats: