Context Navigation

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

Last change on this file since 7682 was 7682, checked in by hudson, 13 years ago
[Ticket 344] - produce asc when there is a hole in the model. Bug fix: holes no longer cause an exception, and are returned as outside points.
File size: 50.2 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 ??
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	inside_poly_indices, outside_poly_indices = \
467	in_and_outside_polygon(point_coordinates,
468	self.mesh.get_boundary_polygon(),
469	closed=True, verbose=verbose)
470
471	# Build n x m interpolation matrix
472	if verbose and len(outside_poly_indices) > 0:
473	log.critical('WARNING: Points outside mesh boundary.')
474
475	# Since you can block, throw a warning, not an error.
476	if verbose and 0 == len(inside_poly_indices):
477	log.critical('WARNING: No points within the mesh!')
478
479	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
480	n = point_coordinates.shape[0] # Nbr of data points
481
482	if verbose: log.critical('Number of datapoints: %d' % n)
483	if verbose: log.critical('Number of basis functions: %d' % m)
484
485	A = Sparse(n,m)
486
487	n = len(inside_poly_indices)
488
489	centroids = []
490
491	# Compute matrix elements for points inside the mesh
492	if verbose: log.critical('Building interpolation matrix from %d points'
493	% n)
494
495	for d, i in enumerate(inside_poly_indices):
496	# For each data_coordinate point
497	if verbose and d%((n+10)/10)==0: log.critical('Doing %d of %d'
498	%(d, n))
499
500	x = point_coordinates[i]
501	element_found, sigma0, sigma1, sigma2, k = \
502	search_tree_of_vertices(self.root, x)
503
504	# Update interpolation matrix A if necessary
505	if element_found is True:
506	# Assign values to matrix A
507	j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0
508	j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
509	j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
510	js = [j0, j1, j2]
511
512	if output_centroids is False:
513	# Weight each vertex according to its distance from x
514	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
515	for j in js:
516	A[i, j] = sigmas[j]
517	else:
518	# If centroids are needed, weight all 3 vertices equally
519	for j in js:
520	A[i, j] = 1.0/3.0
521	centroids.append(self.mesh.centroid_coordinates[k])
522	else:
523	# Mesh has a hole - move this point to the outside list
524	# num.delete(inside_poly_indices,[i], axis=0)
525	num.append(outside_poly_indices, [i], axis=0)
526	# msg = 'Could not find triangle for point', x
527	# raise Exception(msg)
528
529	return A, inside_poly_indices, outside_poly_indices, centroids
530
531
532
533
534
535
536	##
537	# @brief ??
538	# @param vertices ??
539	# @param vertex_attributes ??
540	# @param triangles ??
541	# @param points ??
542	# @param max_points_per_cell ??
543	# @param start_blocking_len ??
544	# @param mesh_origin ??
545	def benchmark_interpolate(vertices,
546	vertex_attributes,
547	triangles, points,
548	max_points_per_cell=None,
549	start_blocking_len=500000,
550	mesh_origin=None):
551	"""
552	points: Interpolate mesh data to these positions.
553	List of coordinate pairs [x, y] of
554	data points or an nx2 numeric array or a Geospatial_data object
555
556	No test for this yet.
557	Note, this has no time the input data has no time dimension. Which is
558	different from most of the data we interpolate, eg sww info.
559
560	Output:
561	Interpolated values at inputted points.
562	"""
563
564	interp = Interpolate(vertices,
565	triangles,
566	max_vertices_per_cell=max_points_per_cell,
567	mesh_origin=mesh_origin)
568
569	calc = interp.interpolate(vertex_attributes,
570	points,
571	start_blocking_len=start_blocking_len)
572
573
574	##
575	# @brief Interpolate quantities at given locations (from .SWW file).
576	# @param sww_file Input .SWW file.
577	# @param points A list of the 'gauges' x,y location.
578	# @param depth_file The name of the output depth file.
579	# @param velocity_x_file Name of the output x velocity file.
580	# @param velocity_y_file Name of the output y velocity file.
581	# @param stage_file Name of the output stage file.
582	# @param froude_file
583	# @param time_thinning Time thinning step to use.
584	# @param verbose True if this function is to be verbose.
585	# @param use_cache True if we are caching.
586	def interpolate_sww2csv(sww_file,
587	points,
588	depth_file,
589	velocity_x_file,
590	velocity_y_file,
591	stage_file=None,
592	froude_file=None,
593	time_thinning=1,
594	verbose=True,
595	use_cache = True):
596	"""
597	Interpolate the quantities at a given set of locations, given
598	an sww file.
599	The results are written to csv files.
600
601	sww_file is the input sww file.
602	points is a list of the 'gauges' x,y location.
603	depth_file is the name of the output depth file
604	velocity_x_file is the name of the output x velocity file.
605	velocity_y_file is the name of the output y velocity file.
606	stage_file is the name of the output stage file.
607
608	In the csv files columns represents the gauges and each row is a
609	time slice.
610
611	Time_thinning_number controls how many timesteps to use. Only
612	timesteps with index%time_thinning_number == 0 will used, or
613	in other words a value of 3, say, will cause the algorithm to
614	use every third time step.
615
616	In the future let points be a points file.
617	And let the user choose the quantities.
618
619	This is currently quite specific.
620	If it is need to be more general, change things.
621	"""
622
623	quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
624	points = ensure_absolute(points)
625	point_count = len(points)
626	callable_sww = file_function(sww_file,
627	quantities=quantities,
628	interpolation_points=points,
629	verbose=verbose,
630	time_thinning=time_thinning,
631	use_cache=use_cache)
632
633	depth_writer = writer(file(depth_file, "wb"))
634	velocity_x_writer = writer(file(velocity_x_file, "wb"))
635	velocity_y_writer = writer(file(velocity_y_file, "wb"))
636	if stage_file is not None:
637	stage_writer = writer(file(stage_file, "wb"))
638	if froude_file is not None:
639	froude_writer = writer(file(froude_file, "wb"))
640
641	# Write heading
642	heading = [str(x[0])+ ':' + str(x[1]) for x in points]
643	heading.insert(0, "time")
644	depth_writer.writerow(heading)
645	velocity_x_writer.writerow(heading)
646	velocity_y_writer.writerow(heading)
647	if stage_file is not None:
648	stage_writer.writerow(heading)
649	if froude_file is not None:
650	froude_writer.writerow(heading)
651
652	for time in callable_sww.get_time():
653	depths = [time]
654	velocity_xs = [time]
655	velocity_ys = [time]
656	if stage_file is not None:
657	stages = [time]
658	if froude_file is not None:
659	froudes = [time]
660	for point_i, point in enumerate(points):
661	quantities = callable_sww(time,point_i)
662
663	w = quantities[0]
664	z = quantities[1]
665	momentum_x = quantities[2]
666	momentum_y = quantities[3]
667	depth = w - z
668
669	if w == NAN or z == NAN or momentum_x == NAN:
670	velocity_x = NAN
671	else:
672	if depth > 1.e-30: # use epsilon
673	velocity_x = momentum_x / depth #Absolute velocity
674	else:
675	velocity_x = 0
676
677	if w == NAN or z == NAN or momentum_y == NAN:
678	velocity_y = NAN
679	else:
680	if depth > 1.e-30: # use epsilon
681	velocity_y = momentum_y / depth #Absolute velocity
682	else:
683	velocity_y = 0
684
685	if depth < 1.e-30: # use epsilon
686	froude = NAN
687	else:
688	froude = sqrt(velocity_xvelocity_x + velocity_yvelocity_y)/ \
689	sqrt(depth * 9.8066) # gravity m/s/s
690
691	depths.append(depth)
692	velocity_xs.append(velocity_x)
693	velocity_ys.append(velocity_y)
694
695	if stage_file is not None:
696	stages.append(w)
697	if froude_file is not None:
698	froudes.append(froude)
699
700	depth_writer.writerow(depths)
701	velocity_x_writer.writerow(velocity_xs)
702	velocity_y_writer.writerow(velocity_ys)
703
704	if stage_file is not None:
705	stage_writer.writerow(stages)
706	if froude_file is not None:
707	froude_writer.writerow(froudes)
708
709
710	##
711	# @brief
712	class Interpolation_function:
713	"""Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
714	which is interpolated from time series defined at vertices of
715	triangular mesh (such as those stored in sww files)
716
717	Let m be the number of vertices, n the number of triangles
718	and p the number of timesteps.
719	Also, let N be the number of interpolation points.
720
721	Mandatory input
722	time: px1 array of monotonously increasing times (float)
723	quantities: Dictionary of arrays or 1 array (float)
724	The arrays must either have dimensions pxm or mx1.
725	The resulting function will be time dependent in
726	the former case while it will be constant with
727	respect to time in the latter case.
728
729	Optional input:
730	quantity_names: List of keys into the quantities dictionary for
731	imposing a particular order on the output vector.
732	vertex_coordinates: mx2 array of coordinates (float)
733	triangles: nx3 array of indices into vertex_coordinates (int)
734	interpolation_points: Nx2 array of coordinates to be interpolated to
735	verbose: Level of reporting
736
737	The quantities returned by the callable object are specified by
738	the list quantities which must contain the names of the
739	quantities to be returned and also reflect the order, e.g. for
740	the shallow water wave equation, on would have
741	quantities = ['stage', 'xmomentum', 'ymomentum']
742
743	The parameter interpolation_points decides at which points interpolated
744	quantities are to be computed whenever object is called.
745	If None, return average value
746
747	FIXME (Ole): Need to allow vertex coordinates and interpolation points to
748	be geospatial data objects
749
750	(FIXME (Ole): This comment should be removed)
751	Time assumed to be relative to starttime
752	All coordinates assume origin of (0,0) - e.g. georeferencing must be
753	taken care of outside this function
754	"""
755
756	##
757	# @brief ??
758	# @param time ??
759	# @param quantities ??
760	# @param quantity_names ??
761	# @param vertex_coordinates ??
762	# @param triangles ??
763	# @param interpolation_points ??
764	# @param time_thinning ??
765	# @param verbose ??
766	# @param gauge_neighbour_id ??
767	def __init__(self,
768	time,
769	quantities,
770	quantity_names=None,
771	vertex_coordinates=None,
772	triangles=None,
773	interpolation_points=None,
774	time_thinning=1,
775	verbose=False,
776	gauge_neighbour_id=None,
777	output_centroids=False):
778	"""Initialise object and build spatial interpolation if required
779
780	Time_thinning_number controls how many timesteps to use. Only timesteps
781	with index%time_thinning_number == 0 will used, or in other words a
782	value of 3, say, will cause the algorithm to use every third time step.
783	"""
784
785	from anuga.config import time_format
786	import types
787
788	if verbose is True:
789	log.critical('Interpolation_function: input checks')
790
791	# Check temporal info
792	time = ensure_numeric(time)
793
794	if not num.alltrue(time[1:] - time[:-1] >= 0):
795	# This message is time consuming to form due to the conversion of
796	msg = 'Time must be a monotonuosly increasing sequence %s' % time
797	raise Exception, msg
798
799	# Check if quantities is a single array only
800	if type(quantities) != types.DictType:
801	quantities = ensure_numeric(quantities)
802	quantity_names = ['Attribute']
803
804	# Make it a dictionary
805	quantities = {quantity_names[0]: quantities}
806
807	# Use keys if no names are specified
808	if quantity_names is None:
809	quantity_names = quantities.keys()
810
811	# Check spatial info
812	if vertex_coordinates is None:
813	self.spatial = False
814	else:
815	# FIXME (Ole): Try ensure_numeric here -
816	# this function knows nothing about georefering.
817	vertex_coordinates = ensure_absolute(vertex_coordinates)
818
819	if triangles is not None:
820	triangles = ensure_numeric(triangles)
821	self.spatial = True
822
823	if verbose is True:
824	log.critical('Interpolation_function: thinning by %d'
825	% time_thinning)
826
827
828	# Thin timesteps if needed
829	# Note array() is used to make the thinned arrays contiguous in memory
830	self.time = num.array(time[::time_thinning])
831	for name in quantity_names:
832	if len(quantities[name].shape) == 2:
833	quantities[name] = num.array(quantities[name][::time_thinning,:])
834
835	if verbose is True:
836	log.critical('Interpolation_function: precomputing')
837
838	# Save for use with statistics
839	self.quantities_range = {}
840	for name in quantity_names:
841	q = quantities[name][:].flatten()
842	self.quantities_range[name] = [min(q), max(q)]
843
844	self.quantity_names = quantity_names
845	self.vertex_coordinates = vertex_coordinates
846	self.interpolation_points = interpolation_points
847
848	self.index = 0 # Initial time index
849	self.precomputed_values = {}
850	self.centroids = []
851
852	# Precomputed spatial interpolation if requested
853	if interpolation_points is not None:
854	#no longer true. sts files have spatial = True but
855	#if self.spatial is False:
856	# raise 'Triangles and vertex_coordinates must be specified'
857	#
858	try:
859	self.interpolation_points = \
860	interpolation_points = ensure_numeric(interpolation_points)
861	except:
862	msg = 'Interpolation points must be an N x 2 numeric array ' \
863	'or a list of points\n'
864	msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...')
865	raise msg
866
867	# Ensure 'mesh_boundary_polygon' is defined
868	mesh_boundary_polygon = None
869
870	if triangles is not None and vertex_coordinates is not None:
871	# Check that all interpolation points fall within
872	# mesh boundary as defined by triangles and vertex_coordinates.
873	from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
874	from anuga.utilities.polygon import outside_polygon
875
876	# Create temporary mesh object from mesh info passed
877	# into this function.
878	mesh = Mesh(vertex_coordinates, triangles)
879	mesh_boundary_polygon = mesh.get_boundary_polygon()
880
881	indices = outside_polygon(interpolation_points,
882	mesh_boundary_polygon)
883
884	# Record result
885	#self.mesh_boundary_polygon = mesh_boundary_polygon
886	self.indices_outside_mesh = indices
887
888	# Report
889	if len(indices) > 0:
890	msg = 'Interpolation points in Interpolation function fall '
891	msg += 'outside specified mesh. Offending points:\n'
892	out_interp_pts = []
893	for i in indices:
894	msg += '%d: %s\n' % (i, interpolation_points[i])
895	out_interp_pts.append(
896	ensure_numeric(interpolation_points[i]))
897
898	if verbose is True:
899	import sys
900	if sys.platform == 'win32':
901	# FIXME (Ole): Why only Windoze?
902	from anuga.utilities.polygon import plot_polygons
903	title = ('Interpolation points fall '
904	'outside specified mesh')
905	plot_polygons([mesh_boundary_polygon,
906	interpolation_points,
907	out_interp_pts],
908	['line', 'point', 'outside'],
909	figname='points_boundary_out',
910	label=title,
911	verbose=verbose)
912
913	# Joaquim Luis suggested this as an Exception, so
914	# that the user can now what the problem is rather than
915	# looking for NaN's. However, NANs are handy as they can
916	# be ignored leaving good points for continued processing.
917	if verbose:
918	log.critical(msg)
919	#raise Exception(msg)
920
921	elif triangles is None and vertex_coordinates is not None: #jj
922	#Dealing with sts file
923	pass
924	else:
925	raise Exception('Sww file function requires both triangles and '
926	'vertex_coordinates. sts file file function '
927	'requires the latter.')
928
929	# Plot boundary and interpolation points,
930	# but only if if 'mesh_boundary_polygon' has data.
931	if verbose is True and mesh_boundary_polygon is not None:
932	import sys
933	if sys.platform == 'win32':
934	from anuga.utilities.polygon import plot_polygons
935	title = ('Interpolation function: '
936	'Polygon and interpolation points')
937	plot_polygons([mesh_boundary_polygon,
938	interpolation_points],
939	['line', 'point'],
940	figname='points_boundary',
941	label=title,
942	verbose=verbose)
943
944	m = len(self.interpolation_points)
945	p = len(self.time)
946
947	for name in quantity_names:
948	self.precomputed_values[name] = num.zeros((p, m), num.float)
949
950	if verbose is True:
951	log.critical('Build interpolator')
952
953
954	# Build interpolator
955	if triangles is not None and vertex_coordinates is not None:
956	if verbose:
957	msg = 'Building interpolation matrix from source mesh '
958	msg += '(%d vertices, %d triangles)' \
959	% (vertex_coordinates.shape[0],
960	triangles.shape[0])
961	log.critical(msg)
962
963	# This one is no longer needed for STS files
964	interpol = Interpolate(vertex_coordinates,
965	triangles,
966	verbose=verbose)
967
968	elif triangles is None and vertex_coordinates is not None:
969	if verbose:
970	log.critical('Interpolation from STS file')
971
972
973
974	if verbose:
975	log.critical('Interpolating (%d interpolation points, %d timesteps).'
976	% (self.interpolation_points.shape[0], self.time.shape[0]))
977
978	if time_thinning > 1:
979	log.critical('Timesteps were thinned by a factor of %d'
980	% time_thinning)
981	else:
982	log.critical()
983
984	for i, t in enumerate(self.time):
985	# Interpolate quantities at this timestep
986	if verbose and i%((p+10)/10) == 0:
987	log.critical(' time step %d of %d' % (i, p))
988
989	for name in quantity_names:
990	if len(quantities[name].shape) == 2:
991	Q = quantities[name][i,:] # Quantities at timestep i
992	else:
993	Q = quantities[name][:] # No time dependency
994
995	if verbose and i%((p+10)/10) == 0:
996	log.critical(' quantity %s, size=%d' % (name, len(Q)))
997
998	# Interpolate
999	if triangles is not None and vertex_coordinates is not None:
1000	result = interpol.interpolate(Q,
1001	point_coordinates=\
1002	self.interpolation_points,
1003	verbose=False,
1004	output_centroids=output_centroids)
1005	self.centroids = interpol.centroids
1006	elif triangles is None and vertex_coordinates is not None:
1007	result = interpolate_polyline(Q,
1008	vertex_coordinates,
1009	gauge_neighbour_id,
1010	interpolation_points=\
1011	self.interpolation_points)
1012
1013	#assert len(result), len(interpolation_points)
1014	self.precomputed_values[name][i, :] = result
1015
1016	# Report
1017	if verbose:
1018	log.critical(self.statistics())
1019	else:
1020	# Store quantitites as is
1021	for name in quantity_names:
1022	self.precomputed_values[name] = quantities[name]
1023
1024	##
1025	# @brief Override object representation method.
1026	def __repr__(self):
1027	# return 'Interpolation function (spatio-temporal)'
1028	return self.statistics()
1029
1030	##
1031	# @brief Evaluate interpolation function
1032	# @param t Model time - must lie within existing timesteps.
1033	# @param point_id Index of one of the preprocessed points.
1034	# @param x ??
1035	# @param y ??
1036	# @return ??
1037	def __call__(self, t, point_id=None, x=None, y=None):
1038	"""Evaluate f(t) or f(t, point_id)
1039
1040	Inputs:
1041	t: time - Model time. Must lie within existing timesteps
1042	point_id: index of one of the preprocessed points.
1043
1044	If spatial info is present and all of point_id
1045	are None an exception is raised
1046
1047	If no spatial info is present, point_id arguments are ignored
1048	making f a function of time only.
1049
1050	FIXME: f(t, x, y) x, y could overrided location, point_id ignored
1051	FIXME: point_id could also be a slice
1052	FIXME: What if x and y are vectors?
1053	FIXME: What about f(x,y) without t?
1054	"""
1055
1056	from math import pi, cos, sin, sqrt
1057
1058	if self.spatial is True:
1059	if point_id is None:
1060	if x is None or y is None:
1061	msg = 'Either point_id or x and y must be specified'
1062	raise Exception(msg)
1063	else:
1064	if self.interpolation_points is None:
1065	msg = 'Interpolation_function must be instantiated ' + \
1066	'with a list of interpolation points before ' + \
1067	'parameter point_id can be used'
1068	raise Exception(msg)
1069
1070	msg = 'Time interval [%.16f:%.16f]' % (self.time[0], self.time[-1])
1071	msg += ' does not match model time: %.16f\n' % t
1072	if t < self.time[0]: raise Modeltime_too_early(msg)
1073	if t > self.time[-1]: raise Modeltime_too_late(msg)
1074
1075	oldindex = self.index #Time index
1076
1077	# Find current time slot
1078	while t > self.time[self.index]: self.index += 1
1079	while t < self.time[self.index]: self.index -= 1
1080
1081	if t == self.time[self.index]:
1082	# Protect against case where t == T[-1] (last time)
1083	# - also works in general when t == T[i]
1084	ratio = 0
1085	else:
1086	# t is now between index and index+1
1087	ratio = ((t - self.time[self.index]) /
1088	(self.time[self.index+1] - self.time[self.index]))
1089
1090	# Compute interpolated values
1091	q = num.zeros(len(self.quantity_names), num.float)
1092	for i, name in enumerate(self.quantity_names):
1093	Q = self.precomputed_values[name]
1094
1095	if self.spatial is False:
1096	# If there is no spatial info
1097	assert len(Q.shape) == 1
1098
1099	Q0 = Q[self.index]
1100	if ratio > 0: Q1 = Q[self.index+1]
1101	else:
1102	if x is not None and y is not None:
1103	# Interpolate to x, y
1104	raise 'x,y interpolation not yet implemented'
1105	else:
1106	# Use precomputed point
1107	Q0 = Q[self.index, point_id]
1108	if ratio > 0:
1109	Q1 = Q[self.index+1, point_id]
1110
1111	# Linear temporal interpolation
1112	if ratio > 0:
1113	if Q0 == NAN and Q1 == NAN:
1114	q[i] = Q0
1115	else:
1116	q[i] = Q0 + ratio*(Q1 - Q0)
1117	else:
1118	q[i] = Q0
1119
1120	# Return vector of interpolated values
1121	# FIXME:
1122	if self.spatial is True:
1123	return q
1124	else:
1125	# Replicate q according to x and y
1126	# This is e.g used for Wind_stress
1127	if x is None or y is None:
1128	return q
1129	else:
1130	try:
1131	N = len(x)
1132	except:
1133	return q
1134	else:
1135	# x is a vector - Create one constant column for each value
1136	N = len(x)
1137	assert len(y) == N, 'x and y must have same length'
1138	res = []
1139	for col in q:
1140	res.append(col*num.ones(N, num.float))
1141
1142	return res
1143
1144	##
1145	# @brief Return model time as a vector of timesteps.
1146	def get_time(self):
1147	"""Return model time as a vector of timesteps
1148	"""
1149	return self.time
1150
1151	##
1152	# @brief Output statistics about interpolation_function.
1153	# @return The statistics string.
1154	def statistics(self):
1155	"""Output statistics about interpolation_function
1156	"""
1157
1158	vertex_coordinates = self.vertex_coordinates
1159	interpolation_points = self.interpolation_points
1160	quantity_names = self.quantity_names
1161	#quantities = self.quantities
1162	precomputed_values = self.precomputed_values
1163
1164	x = vertex_coordinates[:,0]
1165	y = vertex_coordinates[:,1]
1166
1167	str = '------------------------------------------------\n'
1168	str += 'Interpolation_function (spatio-temporal) statistics:\n'
1169	str += ' Extent:\n'
1170	str += ' x in [%f, %f], len(x) == %d\n'\
1171	%(min(x), max(x), len(x))
1172	str += ' y in [%f, %f], len(y) == %d\n'\
1173	%(min(y), max(y), len(y))
1174	str += ' t in [%f, %f], len(t) == %d\n'\
1175	%(min(self.time), max(self.time), len(self.time))
1176	str += ' Quantities:\n'
1177	for name in quantity_names:
1178	minq, maxq = self.quantities_range[name]
1179	str += ' %s in [%f, %f]\n' %(name, minq, maxq)
1180	#q = quantities[name][:].flatten()
1181	#str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
1182
1183	if interpolation_points is not None:
1184	str += ' Interpolation points (xi, eta):'\
1185	' number of points == %d\n' %interpolation_points.shape[0]
1186	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
1187	max(interpolation_points[:,0]))
1188	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
1189	max(interpolation_points[:,1]))
1190	str += ' Interpolated quantities (over all timesteps):\n'
1191
1192	for name in quantity_names:
1193	q = precomputed_values[name][:].flatten()
1194	str += ' %s at interpolation points in [%f, %f]\n'\
1195	%(name, min(q), max(q))
1196	str += '------------------------------------------------\n'
1197
1198	return str
1199
1200
1201	##
1202	# @brief ??
1203	# @param sww_file ??
1204	# @param time ??
1205	# @param interpolation_points ??
1206	# @param quantity_names ??
1207	# @param verbose ??
1208	# @note Obsolete. Use file_function() in utils.
1209	def interpolate_sww(sww_file, time, interpolation_points,
1210	quantity_names=None, verbose=False):
1211	"""
1212	obsolete.
1213	use file_function in utils
1214	"""
1215
1216	#open sww file
1217	x, y, volumes, time, quantities = read_sww(sww_file)
1218	log.critical("x=%s" % str(x))
1219	log.critical("y=%s" % str(y))
1220
1221	log.critical("time=%s" % str(time))
1222	log.critical("quantities=%s" % str(quantities))
1223
1224	#Add the x and y together
1225	vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),
1226	axis=1)
1227
1228	#Will return the quantity values at the specified times and locations
1229	interp = Interpolation_interface(time,
1230	quantities,
1231	quantity_names=quantity_names,
1232	vertex_coordinates=vertex_coordinates,
1233	triangles=volumes,
1234	interpolation_points=interpolation_points,
1235	verbose=verbose)
1236
1237
1238	##
1239	# @brief ??
1240	# @param file_name Name of the .SWW file to read.
1241	def read_sww(file_name):
1242	"""
1243	obsolete - Nothing should be calling this
1244
1245	Read in an sww file.
1246
1247	Input;
1248	file_name - the sww file
1249
1250	Output;
1251	x - Vector of x values
1252	y - Vector of y values
1253	z - Vector of bed elevation
1254	volumes - Array. Each row has 3 values, representing
1255	the vertices that define the volume
1256	time - Vector of the times where there is stage information
1257	stage - array with respect to time and vertices (x,y)
1258	"""
1259
1260	msg = 'Function read_sww in interpolat.py is obsolete'
1261	raise Exception, msg
1262
1263	#FIXME Have this reader as part of data_manager?
1264
1265	from Scientific.IO.NetCDF import NetCDFFile
1266	import tempfile
1267	import sys
1268	import os
1269
1270	#Check contents
1271	#Get NetCDF
1272
1273	# see if the file is there. Throw a QUIET IO error if it isn't
1274	# I don't think I could implement the above
1275
1276	#throws prints to screen if file not present
1277	junk = tempfile.mktemp(".txt")
1278	fd = open(junk,'w')
1279	stdout = sys.stdout
1280	sys.stdout = fd
1281	fid = NetCDFFile(file_name, netcdf_mode_r)
1282	sys.stdout = stdout
1283	fd.close()
1284	#clean up
1285	os.remove(junk)
1286
1287	# Get the variables
1288	x = fid.variables['x'][:]
1289	y = fid.variables['y'][:]
1290	volumes = fid.variables['volumes'][:]
1291	time = fid.variables['time'][:]
1292
1293	keys = fid.variables.keys()
1294	keys.remove("x")
1295	keys.remove("y")
1296	keys.remove("volumes")
1297	keys.remove("time")
1298	#Turn NetCDF objects into numeric arrays
1299	quantities = {}
1300	for name in keys:
1301	quantities[name] = fid.variables[name][:]
1302
1303	fid.close()
1304	return x, y, volumes, time, quantities
1305
1306
1307	#-------------------------------------------------------------
1308	if __name__ == "__main__":
1309	names = ["x","y"]
1310	someiterable = [[1,2],[3,4]]
1311	csvwriter = writer(file("some.csv", "wb"))
1312	csvwriter.writerow(names)
1313	for row in someiterable:
1314	csvwriter.writerow(row)

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

Download in other formats: