Context Navigation

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

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

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

Download in other formats: