Context Navigation

source: anuga_core/source/anuga/fit_interpolate/interpolate.py @ 7707

Last change on this file since 7707 was 7707, checked in by hudson, 14 years ago
New quadtree implementation - unoptimised and no tree balancing. A couple of failing unit tests to fix.
File size: 50.7 KB

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

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

Download in other formats: