Context Navigation

source: branches/numpy/anuga/fit_interpolate/interpolate.py @ 6410

Last change on this file since 6410 was 6410, checked in by rwilson, 15 years ago
numpy changes.
File size: 47.6 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 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: 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]), num.int) #array default#
299	else:
300	z = num.zeros((0,), num.int) #array default#
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][:].flatten()
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	title = ('Interpolation points fall '
871	'outside specified mesh')
872	plot_polygons([mesh_boundary_polygon,
873	interpolation_points,
874	out_interp_pts],
875	['line', 'point', 'outside'],
876	figname='points_boundary_out',
877	label=title,
878	verbose=verbose)
879
880	# Joaquim Luis suggested this as an Exception, so
881	# that the user can now what the problem is rather than
882	# looking for NaN's. However, NANs are handy as they can
883	# be ignored leaving good points for continued processing.
884	if verbose:
885	print msg
886	#raise Exception(msg)
887
888	elif triangles is None and vertex_coordinates is not None: #jj
889	#Dealing with sts file
890	pass
891	else:
892	raise Exception('Sww file function requires both triangles and '
893	'vertex_coordinates. sts file file function '
894	'requires the latter.')
895
896	# Plot boundary and interpolation points
897	if verbose is True:
898	import sys
899	if sys.platform == 'win32':
900	from anuga.utilities.polygon import plot_polygons
901	title = ('Interpolation function: '
902	'Polygon and interpolation points')
903	plot_polygons([mesh_boundary_polygon,
904	interpolation_points],
905	['line', 'point'],
906	figname='points_boundary',
907	label=title,
908	verbose=verbose)
909
910	m = len(self.interpolation_points)
911	p = len(self.time)
912
913	for name in quantity_names:
914	self.precomputed_values[name] = num.zeros((p, m), num.float)
915
916	if verbose is True:
917	print 'Build interpolator'
918
919
920	# Build interpolator
921	if triangles is not None and vertex_coordinates is not None:
922	if verbose:
923	msg = 'Building interpolation matrix from source mesh '
924	msg += '(%d vertices, %d triangles)' \
925	% (vertex_coordinates.shape[0],
926	triangles.shape[0])
927	print msg
928
929	# This one is no longer needed for STS files
930	interpol = Interpolate(vertex_coordinates,
931	triangles,
932	verbose=verbose)
933
934	elif triangles is None and vertex_coordinates is not None:
935	if verbose:
936	msg = 'Interpolation from STS file'
937	print msg
938
939
940
941	if verbose:
942	print 'Interpolating (%d interpolation points, %d timesteps).' \
943	%(self.interpolation_points.shape[0], self.time.shape[0]),
944
945	if time_thinning > 1:
946	print 'Timesteps were thinned by a factor of %d' \
947	% time_thinning
948	else:
949	print
950
951	for i, t in enumerate(self.time):
952	# Interpolate quantities at this timestep
953	if verbose and i%((p+10)/10) == 0:
954	print ' time step %d of %d' %(i, p)
955
956	for name in quantity_names:
957	if len(quantities[name].shape) == 2:
958	Q = quantities[name][i,:] # Quantities at timestep i
959	else:
960	Q = quantities[name][:] # No time dependency
961
962	if verbose and i%((p+10)/10) == 0:
963	print ' quantity %s, size=%d' %(name, len(Q))
964
965	# Interpolate
966	if triangles is not None and vertex_coordinates is not None:
967	result = interpol.interpolate(Q,
968	point_coordinates=\
969	self.interpolation_points,
970	verbose=False) # No clutter
971	elif triangles is None and vertex_coordinates is not None:
972	result = interpolate_polyline(Q,
973	vertex_coordinates,
974	gauge_neighbour_id,
975	interpolation_points=\
976	self.interpolation_points)
977
978	#assert len(result), len(interpolation_points)
979	self.precomputed_values[name][i, :] = result
980
981	# Report
982	if verbose:
983	print self.statistics()
984	else:
985	# Store quantitites as is
986	for name in quantity_names:
987	self.precomputed_values[name] = quantities[name]
988
989	##
990	# @brief Override object representation method.
991	def __repr__(self):
992	# return 'Interpolation function (spatio-temporal)'
993	return self.statistics()
994
995	##
996	# @brief Override object() method.
997	# @param t Model time - must lie within existing timesteps.
998	# @param point_id Index of one of the preprocessed points.
999	# @param x ??
1000	# @param y ??
1001	# @return ??
1002	def __call__(self, t, point_id=None, x=None, y=None):
1003	"""Evaluate f(t) or f(t, point_id)
1004
1005	Inputs:
1006	t: time - Model time. Must lie within existing timesteps
1007	point_id: index of one of the preprocessed points.
1008
1009	If spatial info is present and all of point_id
1010	are None an exception is raised
1011
1012	If no spatial info is present, point_id arguments are ignored
1013	making f a function of time only.
1014
1015	FIXME: f(t, x, y) x, y could overrided location, point_id ignored
1016	FIXME: point_id could also be a slice
1017	FIXME: What if x and y are vectors?
1018	FIXME: What about f(x,y) without t?
1019	"""
1020
1021	from math import pi, cos, sin, sqrt
1022	from anuga.abstract_2d_finite_volumes.util import mean
1023
1024	if self.spatial is True:
1025	if point_id is None:
1026	if x is None or y is None:
1027	msg = 'Either point_id or x and y must be specified'
1028	raise Exception(msg)
1029	else:
1030	if self.interpolation_points is None:
1031	msg = 'Interpolation_function must be instantiated ' + \
1032	'with a list of interpolation points before ' + \
1033	'parameter point_id can be used'
1034	raise Exception(msg)
1035
1036	msg = 'Time interval [%.16f:%.16f]' % (self.time[0], self.time[-1])
1037	msg += ' does not match model time: %.16f\n' % t
1038	if t < self.time[0]: raise Modeltime_too_early(msg)
1039	if t > self.time[-1]: raise Modeltime_too_late(msg)
1040
1041	oldindex = self.index #Time index
1042
1043	# Find current time slot
1044	while t > self.time[self.index]: self.index += 1
1045	while t < self.time[self.index]: self.index -= 1
1046
1047	if t == self.time[self.index]:
1048	# Protect against case where t == T[-1] (last time)
1049	# - also works in general when t == T[i]
1050	ratio = 0
1051	else:
1052	# t is now between index and index+1
1053	ratio = ((t - self.time[self.index]) /
1054	(self.time[self.index+1] - self.time[self.index]))
1055
1056	# Compute interpolated values
1057	q = num.zeros(len(self.quantity_names), num.float)
1058	for i, name in enumerate(self.quantity_names):
1059	Q = self.precomputed_values[name]
1060
1061	if self.spatial is False:
1062	# If there is no spatial info
1063	assert len(Q.shape) == 1
1064
1065	Q0 = Q[self.index]
1066	if ratio > 0: Q1 = Q[self.index+1]
1067	else:
1068	if x is not None and y is not None:
1069	# Interpolate to x, y
1070	raise 'x,y interpolation not yet implemented'
1071	else:
1072	# Use precomputed point
1073	Q0 = Q[self.index, point_id]
1074	if ratio > 0:
1075	Q1 = Q[self.index+1, point_id]
1076
1077	# Linear temporal interpolation
1078	if ratio > 0:
1079	if Q0 == NAN and Q1 == NAN:
1080	q[i] = Q0
1081	else:
1082	q[i] = Q0 + ratio*(Q1 - Q0)
1083	else:
1084	q[i] = Q0
1085
1086	# Return vector of interpolated values
1087	# FIXME:
1088	if self.spatial is True:
1089	return q
1090	else:
1091	# Replicate q according to x and y
1092	# This is e.g used for Wind_stress
1093	if x is None or y is None:
1094	return q
1095	else:
1096	try:
1097	N = len(x)
1098	except:
1099	return q
1100	else:
1101	# x is a vector - Create one constant column for each value
1102	N = len(x)
1103	assert len(y) == N, 'x and y must have same length'
1104	res = []
1105	for col in q:
1106	res.append(col*num.ones(N, num.float))
1107
1108	return res
1109
1110	##
1111	# @brief Return model time as a vector of timesteps.
1112	def get_time(self):
1113	"""Return model time as a vector of timesteps
1114	"""
1115	return self.time
1116
1117	##
1118	# @brief Output statistics about interpolation_function.
1119	# @return The statistics string.
1120	def statistics(self):
1121	"""Output statistics about interpolation_function
1122	"""
1123
1124	vertex_coordinates = self.vertex_coordinates
1125	interpolation_points = self.interpolation_points
1126	quantity_names = self.quantity_names
1127	#quantities = self.quantities
1128	precomputed_values = self.precomputed_values
1129
1130	x = vertex_coordinates[:,0]
1131	y = vertex_coordinates[:,1]
1132
1133	str = '------------------------------------------------\n'
1134	str += 'Interpolation_function (spatio-temporal) statistics:\n'
1135	str += ' Extent:\n'
1136	str += ' x in [%f, %f], len(x) == %d\n'\
1137	%(min(x), max(x), len(x))
1138	str += ' y in [%f, %f], len(y) == %d\n'\
1139	%(min(y), max(y), len(y))
1140	str += ' t in [%f, %f], len(t) == %d\n'\
1141	%(min(self.time), max(self.time), len(self.time))
1142	str += ' Quantities:\n'
1143	for name in quantity_names:
1144	minq, maxq = self.quantities_range[name]
1145	str += ' %s in [%f, %f]\n' %(name, minq, maxq)
1146	#q = quantities[name][:].flatten()
1147	#str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
1148
1149	if interpolation_points is not None:
1150	str += ' Interpolation points (xi, eta):'\
1151	' number of points == %d\n' %interpolation_points.shape[0]
1152	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
1153	max(interpolation_points[:,0]))
1154	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
1155	max(interpolation_points[:,1]))
1156	str += ' Interpolated quantities (over all timesteps):\n'
1157
1158	for name in quantity_names:
1159	q = precomputed_values[name][:].flatten()
1160	str += ' %s at interpolation points in [%f, %f]\n'\
1161	%(name, min(q), max(q))
1162	str += '------------------------------------------------\n'
1163
1164	return str
1165
1166
1167	##
1168	# @brief ??
1169	# @param sww_file ??
1170	# @param time ??
1171	# @param interpolation_points ??
1172	# @param quantity_names ??
1173	# @param verbose ??
1174	# @note Obsolete. Use file_function() in utils.
1175	def interpolate_sww(sww_file, time, interpolation_points,
1176	quantity_names=None, verbose=False):
1177	"""
1178	obsolete.
1179	use file_function in utils
1180	"""
1181
1182	#open sww file
1183	x, y, volumes, time, quantities = read_sww(sww_file)
1184	print "x",x
1185	print "y",y
1186
1187	print "time", time
1188	print "quantities", quantities
1189
1190	#Add the x and y together
1191	vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),axis=1)
1192
1193	#Will return the quantity values at the specified times and locations
1194	interp = Interpolation_interface(time,
1195	quantities,
1196	quantity_names=quantity_names,
1197	vertex_coordinates=vertex_coordinates,
1198	triangles=volumes,
1199	interpolation_points=interpolation_points,
1200	verbose=verbose)
1201
1202
1203	##
1204	# @brief ??
1205	# @param file_name Name of the .SWW file to read.
1206	def read_sww(file_name):
1207	"""
1208	obsolete - Nothing should be calling this
1209
1210	Read in an sww file.
1211
1212	Input;
1213	file_name - the sww file
1214
1215	Output;
1216	x - Vector of x values
1217	y - Vector of y values
1218	z - Vector of bed elevation
1219	volumes - Array. Each row has 3 values, representing
1220	the vertices that define the volume
1221	time - Vector of the times where there is stage information
1222	stage - array with respect to time and vertices (x,y)
1223	"""
1224
1225	msg = 'Function read_sww in interpolat.py is obsolete'
1226	raise Exception, msg
1227
1228	#FIXME Have this reader as part of data_manager?
1229
1230	from Scientific.IO.NetCDF import NetCDFFile
1231	import tempfile
1232	import sys
1233	import os
1234
1235	#Check contents
1236	#Get NetCDF
1237
1238	# see if the file is there. Throw a QUIET IO error if it isn't
1239	# I don't think I could implement the above
1240
1241	#throws prints to screen if file not present
1242	junk = tempfile.mktemp(".txt")
1243	fd = open(junk,'w')
1244	stdout = sys.stdout
1245	sys.stdout = fd
1246	fid = NetCDFFile(file_name, netcdf_mode_r)
1247	sys.stdout = stdout
1248	fd.close()
1249	#clean up
1250	os.remove(junk)
1251
1252	# Get the variables
1253	x = fid.variables['x'][:]
1254	y = fid.variables['y'][:]
1255	volumes = fid.variables['volumes'][:]
1256	time = fid.variables['time'][:]
1257
1258	keys = fid.variables.keys()
1259	keys.remove("x")
1260	keys.remove("y")
1261	keys.remove("volumes")
1262	keys.remove("time")
1263	#Turn NetCDF objects into numeric arrays
1264	quantities = {}
1265	for name in keys:
1266	quantities[name] = fid.variables[name][:]
1267
1268	fid.close()
1269	return x, y, volumes, time, quantities
1270
1271
1272	#-------------------------------------------------------------
1273	if __name__ == "__main__":
1274	names = ["x","y"]
1275	someiterable = [[1,2],[3,4]]
1276	csvwriter = writer(file("some.csv", "wb"))
1277	csvwriter.writerow(names)
1278	for row in someiterable:
1279	csvwriter.writerow(row)

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

Download in other formats: