Context Navigation

source: branches/Numeric_anuga_source/anuga/fit_interpolate/interpolate.py @ 7258

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

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

Download in other formats: