Context Navigation

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

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

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

Download in other formats: