Context Navigation

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

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

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

Download in other formats: