Context Navigation

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

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

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

Download in other formats: