Context Navigation

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

Last change on this file since 4820 was 4820, checked in by ole, 16 years ago
Addressed ticket:217
File size: 33.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	from warnings import warn
23	from math import sqrt
24	from csv import writer, DictWriter
25
26	from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
27	ArrayType, allclose, take, NewAxis, arange
28
29	from anuga.caching.caching import cache
30	from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
31	from anuga.utilities.sparse import Sparse, Sparse_CSR
32	from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError
33	from anuga.coordinate_transforms.geo_reference import Geo_reference
34	from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
35	from anuga.utilities.polygon import in_and_outside_polygon
36	from anuga.geospatial_data.geospatial_data import Geospatial_data
37	from anuga.geospatial_data.geospatial_data import ensure_absolute
38	from anuga.fit_interpolate.search_functions import search_tree_of_vertices
39	from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
40	from anuga.abstract_2d_finite_volumes.util import file_function
41
42
43	class Interpolate (FitInterpolate):
44
45	def __init__(self,
46	vertex_coordinates,
47	triangles,
48	mesh_origin=None,
49	verbose=False,
50	max_vertices_per_cell=None):
51
52
53	""" Build interpolation matrix mapping from
54	function values at vertices to function values at data points
55
56	Inputs:
57
58	vertex_coordinates: List of coordinate pairs [xi, eta] of
59	points constituting a mesh (or an m x 2 Numeric array or
60	a geospatial object)
61	Points may appear multiple times
62	(e.g. if vertices have discontinuities)
63
64	triangles: List of 3-tuples (or a Numeric array) of
65	integers representing indices of all vertices in the mesh.
66
67	mesh_origin: A geo_reference object or 3-tuples consisting of
68	UTM zone, easting and northing.
69	If specified vertex coordinates are assumed to be
70	relative to their respective origins.
71
72	max_vertices_per_cell: Number of vertices in a quad tree cell
73	at which the cell is split into 4.
74
75	Note: Don't supply a vertex coords as a geospatial object and
76	a mesh origin, since geospatial has its own mesh origin.
77	"""
78
79	# FIXME (Ole): Need an input check
80
81	# Initialise variabels
82	self._A_can_be_reused = False
83	self._point_coordinates = None
84
85	FitInterpolate.__init__(self,
86	vertex_coordinates=vertex_coordinates,
87	triangles=triangles,
88	mesh_origin=mesh_origin,
89	verbose=verbose,
90	max_vertices_per_cell=max_vertices_per_cell)
91
92	# FIXME: What is a good start_blocking_len value?
93	def interpolate(self,
94	f,
95	point_coordinates=None,
96	start_blocking_len=500000,
97	verbose=False):
98	"""Interpolate mesh data f to determine values, z, at points.
99
100	f is the data on the mesh vertices.
101
102	The mesh values representing a smooth surface are
103	assumed to be specified in f.
104
105	Inputs:
106	f: Vector or array of data at the mesh vertices.
107	If f is an array, interpolation will be done for each column as
108	per underlying matrix-matrix multiplication
109	point_coordinates: Interpolate mesh data to these positions.
110	List of coordinate pairs [x, y] of
111	data points or an nx2 Numeric array or a Geospatial_data object
112
113	If point_coordinates is absent, the points inputted last time
114	this method was called are used, if possible.
115	start_blocking_len: If the # of points is more or greater than this,
116	start blocking
117
118	Output:
119	Interpolated values at inputted points (z).
120	"""
121
122
123	# FIXME (Ole): Need an input check that dimensions are compatible
124
125	# FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
126	# method is called even if interpolation points are unchanged.
127
128	#print "point_coordinates interpolate.interpolate", point_coordinates
129	if isinstance(point_coordinates, Geospatial_data):
130	point_coordinates = point_coordinates.get_data_points( \
131	absolute = True)
132
133	# Can I interpolate, based on previous point_coordinates?
134	if point_coordinates is None:
135	if self._A_can_be_reused is True and \
136	len(self._point_coordinates) < start_blocking_len:
137	z = self._get_point_data_z(f,
138	verbose=verbose)
139	elif self._point_coordinates is not None:
140	# if verbose, give warning
141	if verbose:
142	print 'WARNING: Recalculating A matrix, due to blocking.'
143	point_coordinates = self._point_coordinates
144	else:
145	#There are no good point_coordinates. import sys; sys.exit()
146	msg = 'ERROR (interpolate.py): No point_coordinates inputted'
147	raise Exception(msg)
148
149	if point_coordinates is not None:
150	self._point_coordinates = point_coordinates
151	if len(point_coordinates) < start_blocking_len or \
152	start_blocking_len == 0:
153	self._A_can_be_reused = True
154	z = self.interpolate_block(f, point_coordinates,
155	verbose=verbose)
156	else:
157	#print 'BLOCKING'
158	#Handle blocking
159	self._A_can_be_reused = False
160	start = 0
161	# creating a dummy array to concatenate to.
162
163	f = ensure_numeric(f, Float)
164	#print "f.shape",f.shape
165	if len(f.shape) > 1:
166	z = zeros((0,f.shape[1]))
167	else:
168	z = zeros((0,))
169
170	for end in range(start_blocking_len,
171	len(point_coordinates),
172	start_blocking_len):
173
174	t = self.interpolate_block(f, point_coordinates[start:end],
175	verbose=verbose)
176	#print "t", t
177	#print "z", z
178	z = concatenate((z,t))
179	start = end
180
181	end = len(point_coordinates)
182	t = self.interpolate_block(f, point_coordinates[start:end],
183	verbose=verbose)
184	z = concatenate((z,t))
185	return z
186
187	def interpolate_block(self, f, point_coordinates=None, verbose=False):
188	"""
189	Call this if you want to control the blocking or make sure blocking
190	doesn't occur.
191
192	Return the point data, z.
193
194	See interpolate for doc info.
195	"""
196	if isinstance(point_coordinates,Geospatial_data):
197	point_coordinates = point_coordinates.get_data_points( \
198	absolute = True)
199	if point_coordinates is not None:
200	self._A = self._build_interpolation_matrix_A(point_coordinates,
201	verbose=verbose)
202	return self._get_point_data_z(f)
203
204	def _get_point_data_z(self, f, verbose=False):
205	"""
206	Return the point data, z.
207
208	Precondition,
209	The _A matrix has been created
210	"""
211	z = self._A * f
212	# Taking into account points outside the mesh.
213	#print "self.outside_poly_indices", self.outside_poly_indices
214	#print "self.inside_poly_indices", self.inside_poly_indices
215	#print "z", z
216	for i in self.outside_poly_indices:
217	z[i] = NAN
218	return z
219
220	def _build_interpolation_matrix_A(self,
221	point_coordinates,
222	verbose=False):
223	"""Build n x m interpolation matrix, where
224	n is the number of data points and
225	m is the number of basis functions phi_k (one per vertex)
226
227	This algorithm uses a quad tree data structure for fast binning
228	of data points
229	origin is a 3-tuple consisting of UTM zone, easting and northing.
230	If specified coordinates are assumed to be relative to this origin.
231
232	This one will override any data_origin that may be specified in
233	instance interpolation
234
235	Preconditions
236	Point_coordindates and mesh vertices have the same origin.
237	"""
238
239	#print 'Building interpolation matrix'
240
241	#Convert point_coordinates to Numeric arrays, in case it was a list.
242	point_coordinates = ensure_numeric(point_coordinates, Float)
243
244	if verbose: print 'Getting indices inside mesh boundary'
245	self.inside_poly_indices, self.outside_poly_indices = \
246	in_and_outside_polygon(point_coordinates,
247	self.mesh.get_boundary_polygon(),
248	closed = True, verbose = verbose)
249
250	#Build n x m interpolation matrix
251	if verbose and len(self.outside_poly_indices) > 0:
252	print '\n WARNING: Points outside mesh boundary. \n'
253	# Since you can block, throw a warning, not an error.
254	if verbose and 0 == len(self.inside_poly_indices):
255	print '\n WARNING: No points within the mesh! \n'
256
257	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
258	n = point_coordinates.shape[0] # Nbr of data points
259
260	if verbose: print 'Number of datapoints: %d' %n
261	if verbose: print 'Number of basis functions: %d' %m
262
263	A = Sparse(n,m)
264
265	n = len(self.inside_poly_indices)
266	#Compute matrix elements for points inside the mesh
267	if verbose: print 'Building interpolation matrix from %d points' %n
268	for d, i in enumerate(self.inside_poly_indices):
269	# For each data_coordinate point
270	if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
271	x = point_coordinates[i]
272	element_found, sigma0, sigma1, sigma2, k = \
273	search_tree_of_vertices(self.root, self.mesh, x)
274
275	# Update interpolation matrix A if necessary
276	if element_found is True:
277	# Assign values to matrix A
278
279	j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0
280	j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
281	j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
282
283	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
284	js = [j0,j1,j2]
285
286	for j in js:
287	A[i,j] = sigmas[j]
288	else:
289	msg = 'Could not find triangle for point', x
290	raise Exception(msg)
291	return A
292
293	def benchmark_interpolate(vertices,
294	vertex_attributes,
295	triangles, points,
296	max_points_per_cell=None,
297	start_blocking_len=500000):
298	"""
299	No test for this yet.
300	Note, this has no time the input data has no time dimension. Which is
301	different from most of the data we interpolate, eg sww info.
302
303	Output:
304	Interpolated values at inputted points.
305	"""
306	interp = Interpolate(vertices,
307	triangles,
308	max_vertices_per_cell=max_points_per_cell)
309
310	calc = interp.interpolate(vertex_attributes
311	,points
312	,start_blocking_len=start_blocking_len)
313	def interpolate_sww2csv(sww_file,
314	points,
315	depth_file,
316	velocity_x_file,
317	velocity_y_file,
318	#quantities = ['depth', 'velocity'],
319	verbose=True,
320	use_cache = True):
321	"""
322	Interpolate the quantities at a given set of locations, given
323	an sww file.
324	The results are written to a csv file.
325
326	In the future let points be a points file.
327	And the user choose the quantities.
328
329	This is currently quite specific.
330	If it need to be more general, chagne things.
331
332	This is really returning speed, not velocity.
333	"""
334	quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
335	#print "points",points
336	points = ensure_absolute(points)
337	point_count = len(points)
338	callable_sww = file_function(sww_file,
339	quantities=quantities,
340	interpolation_points=points,
341	verbose=verbose,
342	use_cache=use_cache)
343
344	depth_writer = writer(file(depth_file, "wb"))
345	velocity_x_writer = writer(file(velocity_x_file, "wb"))
346	velocity_y_writer = writer(file(velocity_y_file, "wb"))
347	# Write heading
348	heading = [str(x[0])+ ':' + str(x[1]) for x in points]
349	heading.insert(0, "time")
350	depth_writer.writerow(heading)
351	velocity_x_writer.writerow(heading)
352	velocity_y_writer.writerow(heading)
353
354	for time in callable_sww.get_time():
355	depths = [time]
356	velocity_xs = [time]
357	velocity_ys = [time]
358	for point_i, point in enumerate(points):
359	quantities = callable_sww(time,point_i)
360	#print "quantities", quantities
361
362	w = quantities[0]
363	z = quantities[1]
364	momentum_x = quantities[2]
365	momentum_y = quantities[3]
366	depth = w - z
367
368	if w == NAN or z == NAN or momentum_x == NAN:
369	velocity_x = NAN
370	else:
371	if depth > 1.e-30: # use epsilon
372	velocity_x = momentum_x / depth #Absolute velocity
373	else:
374	velocity_x = 0
375	if w == NAN or z == NAN or momentum_y == NAN:
376	velocity_y = NAN
377	else:
378	if depth > 1.e-30: # use epsilon
379	velocity_y = momentum_y / depth #Absolute velocity
380	else:
381	velocity_y = 0
382	depths.append(depth)
383	velocity_xs.append(velocity_x)
384	velocity_ys.append(velocity_y)
385	depth_writer.writerow(depths)
386	velocity_x_writer.writerow(velocity_xs)
387	velocity_y_writer.writerow(velocity_ys)
388
389
390	class Interpolation_function:
391	"""Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
392	which is interpolated from time series defined at vertices of
393	triangular mesh (such as those stored in sww files)
394
395	Let m be the number of vertices, n the number of triangles
396	and p the number of timesteps.
397
398	Mandatory input
399	time: px1 array of monotonously increasing times (Float)
400	quantities: Dictionary of arrays or 1 array (Float)
401	The arrays must either have dimensions pxm or mx1.
402	The resulting function will be time dependent in
403	the former case while it will be constan with
404	respect to time in the latter case.
405
406	Optional input:
407	quantity_names: List of keys into the quantities dictionary
408	vertex_coordinates: mx2 array of coordinates (Float)
409	triangles: nx3 array of indices into vertex_coordinates (Int)
410	interpolation_points: Nx2 array of coordinates to be interpolated to
411	verbose: Level of reporting
412
413
414	The quantities returned by the callable object are specified by
415	the list quantities which must contain the names of the
416	quantities to be returned and also reflect the order, e.g. for
417	the shallow water wave equation, on would have
418	quantities = ['stage', 'xmomentum', 'ymomentum']
419
420	The parameter interpolation_points decides at which points interpolated
421	quantities are to be computed whenever object is called.
422	If None, return average value
423
424	FIXME (Ole): Need to allow vertex coordinates and interpolation points to be
425	geospatial data objects
426
427	Time assumed to be relative to starttime
428	"""
429
430
431	def __init__(self,
432	time,
433	quantities,
434	quantity_names=None,
435	vertex_coordinates=None,
436	triangles=None,
437	interpolation_points=None,
438	time_thinning=1,
439	verbose=False):
440	"""Initialise object and build spatial interpolation if required
441
442	Time_thinning_number controls how many timesteps to use. Only timesteps with
443	index%time_thinning_number == 0 will used, or in other words a value of 3, say,
444	will cause the algorithm to use every third time step.
445	"""
446
447	from Numeric import array, zeros, Float, alltrue, concatenate,\
448	reshape, ArrayType
449
450
451	from anuga.config import time_format
452	import types
453
454
455	# Check temporal info
456	time = ensure_numeric(time)
457	msg = 'Time must be a monotonuosly '
458	msg += 'increasing sequence %s' %time
459	assert alltrue(time[1:] - time[:-1] >= 0 ), msg
460
461
462	# Check if quantities is a single array only
463	if type(quantities) != types.DictType:
464	quantities = ensure_numeric(quantities)
465	quantity_names = ['Attribute']
466
467	#Make it a dictionary
468	quantities = {quantity_names[0]: quantities}
469
470
471	# Use keys if no names are specified
472	if quantity_names is None:
473	quantity_names = quantities.keys()
474
475
476	# Check spatial info
477	if vertex_coordinates is None:
478	self.spatial = False
479	else:
480	vertex_coordinates = ensure_absolute(vertex_coordinates)
481
482	assert triangles is not None, 'Triangles array must be specified'
483	triangles = ensure_numeric(triangles)
484	self.spatial = True
485
486	# Thin timesteps if needed
487	# Note array() is used to make the thinned arrays contiguous in memory
488	self.time = array(time[::time_thinning])
489	for name in quantity_names:
490	if len(quantities[name].shape) == 2:
491	quantities[name] = array(quantities[name][::time_thinning,:])
492
493	# Save for use with statistics
494	self.quantities_range = {}
495	for name in quantity_names:
496	q = quantities[name][:].flat
497	self.quantities_range[name] = [min(q), max(q)]
498
499	self.quantity_names = quantity_names
500	self.vertex_coordinates = vertex_coordinates
501	self.interpolation_points = interpolation_points
502
503
504	self.index = 0 # Initial time index
505	self.precomputed_values = {}
506
507
508	# Precomputed spatial interpolation if requested
509	if interpolation_points is not None:
510	if self.spatial is False:
511	raise 'Triangles and vertex_coordinates must be specified'
512
513	try:
514	self.interpolation_points = ensure_numeric(interpolation_points)
515	except:
516	msg = 'Interpolation points must be an N x 2 Numeric array '+\
517	'or a list of points\n'
518	msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
519	'...')
520	raise msg
521
522
523	# Check that all interpolation points fall within
524	# mesh boundary as defined by triangles and vertex_coordinates.
525	from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
526	from anuga.utilities.polygon import outside_polygon
527
528	mesh = Mesh(vertex_coordinates, triangles)
529
530	indices = outside_polygon(interpolation_points,
531	mesh.get_boundary_polygon())
532
533
534	# Record result
535	self.indices_outside_mesh = indices
536
537	# Report
538	if len(indices) > 0:
539	msg = 'Interpolation points in Interpolation function fall '
540	msg += 'outside specified mesh. '
541	msg += 'Offending points:\n'
542	for i in indices:
543	msg += '%d: %s\n' %(i, interpolation_points[i])
544
545	# Joaquim Luis suggested this as an Exception, so
546	# that the user can now what the problem is rather than
547	# looking for NaN's. However, NANs are handy as they can
548	# be ignored leaving good points for continued processing.
549	if verbose:
550	print msg
551	#raise Exception(msg)
552
553
554
555
556	m = len(self.interpolation_points)
557	p = len(self.time)
558
559	for name in quantity_names:
560	self.precomputed_values[name] = zeros((p, m), Float)
561
562	# Build interpolator
563	if verbose:
564	msg = 'Building interpolation matrix from source mesh '
565	msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
566	triangles.shape[0])
567	print msg
568
569
570	interpol = Interpolate(vertex_coordinates,
571	triangles,
572	verbose=verbose)
573
574	if verbose:
575	print 'Interpolating (%d interpolation points, %d timesteps).'\
576	%(self.interpolation_points.shape[0], self.time.shape[0]),
577
578	if time_thinning > 1:
579	print 'Timesteps were thinned by a factor of %d' %time_thinning
580	else:
581	print
582
583	for i, t in enumerate(self.time):
584	# Interpolate quantities at this timestep
585	if verbose and i%((p+10)/10)==0:
586	print ' time step %d of %d' %(i, p)
587
588	for name in quantity_names:
589	if len(quantities[name].shape) == 2:
590	Q = quantities[name][i,:] # Quantities at timestep i
591	else:
592	Q = quantities[name][:] # No time dependency
593
594	# Interpolate
595	result = interpol.interpolate(Q,
596	point_coordinates=\
597	self.interpolation_points)
598
599	#assert len(result), len(interpolation_points)
600	self.precomputed_values[name][i, :] = result
601
602
603	# Report
604	if verbose:
605	print self.statistics()
606	#self.print_statistics()
607
608	else:
609	# Store quantitites as is
610	for name in quantity_names:
611	self.precomputed_values[name] = quantities[name]
612
613	def __repr__(self):
614	# return 'Interpolation function (spatio-temporal)'
615	return self.statistics()
616
617	def __call__(self, t, point_id=None, x=None, y=None):
618	"""Evaluate f(t) or f(t, point_id)
619
620	Inputs:
621	t: time - Model time. Must lie within existing timesteps
622	point_id: index of one of the preprocessed points.
623
624
625	If spatial info is present and all of point_id
626	are None an exception is raised
627
628	If no spatial info is present, point_id arguments are ignored
629	making f a function of time only.
630
631
632	FIXME: f(t, x, y) x, y could overrided location, point_id ignored
633	FIXME: point_id could also be a slice
634	FIXME: What if x and y are vectors?
635	FIXME: What about f(x,y) without t?
636	"""
637
638	from math import pi, cos, sin, sqrt
639	from Numeric import zeros, Float
640	from anuga.abstract_2d_finite_volumes.util import mean
641
642	if self.spatial is True:
643	if point_id is None:
644	if x is None or y is None:
645	msg = 'Either point_id or x and y must be specified'
646	raise Exception(msg)
647	else:
648	if self.interpolation_points is None:
649	msg = 'Interpolation_function must be instantiated ' +\
650	'with a list of interpolation points before parameter ' +\
651	'point_id can be used'
652	raise Exception(msg)
653
654	msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
655	msg += ' does not match model time: %.16f\n' %t
656	if t < self.time[0]: raise Exception(msg)
657	if t > self.time[-1]: raise Exception(msg)
658
659	oldindex = self.index #Time index
660
661	# Find current time slot
662	while t > self.time[self.index]: self.index += 1
663	while t < self.time[self.index]: self.index -= 1
664
665	if t == self.time[self.index]:
666	# Protect against case where t == T[-1] (last time)
667	# - also works in general when t == T[i]
668	ratio = 0
669	else:
670	# t is now between index and index+1
671	ratio = (t - self.time[self.index])/\
672	(self.time[self.index+1] - self.time[self.index])
673
674	# Compute interpolated values
675	q = zeros(len(self.quantity_names), Float)
676	# print "self.precomputed_values", self.precomputed_values
677	for i, name in enumerate(self.quantity_names):
678	Q = self.precomputed_values[name]
679
680	if self.spatial is False:
681	# If there is no spatial info
682	assert len(Q.shape) == 1
683
684	Q0 = Q[self.index]
685	if ratio > 0: Q1 = Q[self.index+1]
686
687	else:
688	if x is not None and y is not None:
689	# Interpolate to x, y
690
691	raise 'x,y interpolation not yet implemented'
692	else:
693	# Use precomputed point
694	Q0 = Q[self.index, point_id]
695	if ratio > 0:
696	Q1 = Q[self.index+1, point_id]
697
698	# Linear temporal interpolation
699	if ratio > 0:
700	if Q0 == NAN and Q1 == NAN:
701	q[i] = Q0
702	else:
703	q[i] = Q0 + ratio*(Q1 - Q0)
704	else:
705	q[i] = Q0
706
707
708	# Return vector of interpolated values
709	# if len(q) == 1:
710	# return q[0]
711	# else:
712	# return q
713
714
715	# Return vector of interpolated values
716	# FIXME:
717	if self.spatial is True:
718	return q
719	else:
720	# Replicate q according to x and y
721	# This is e.g used for Wind_stress
722	if x is None or y is None:
723	return q
724	else:
725	try:
726	N = len(x)
727	except:
728	return q
729	else:
730	from Numeric import ones, Float
731	# x is a vector - Create one constant column for each value
732	N = len(x)
733	assert len(y) == N, 'x and y must have same length'
734	res = []
735	for col in q:
736	res.append(col*ones(N, Float))
737
738	return res
739
740
741	def get_time(self):
742	"""Return model time as a vector of timesteps
743	"""
744	return self.time
745
746
747	def statistics(self):
748	"""Output statistics about interpolation_function
749	"""
750
751	vertex_coordinates = self.vertex_coordinates
752	interpolation_points = self.interpolation_points
753	quantity_names = self.quantity_names
754	#quantities = self.quantities
755	precomputed_values = self.precomputed_values
756
757	x = vertex_coordinates[:,0]
758	y = vertex_coordinates[:,1]
759
760	str = '------------------------------------------------\n'
761	str += 'Interpolation_function (spatio-temporal) statistics:\n'
762	str += ' Extent:\n'
763	str += ' x in [%f, %f], len(x) == %d\n'\
764	%(min(x), max(x), len(x))
765	str += ' y in [%f, %f], len(y) == %d\n'\
766	%(min(y), max(y), len(y))
767	str += ' t in [%f, %f], len(t) == %d\n'\
768	%(min(self.time), max(self.time), len(self.time))
769	str += ' Quantities:\n'
770	for name in quantity_names:
771	minq, maxq = self.quantities_range[name]
772	str += ' %s in [%f, %f]\n' %(name, minq, maxq)
773	#q = quantities[name][:].flat
774	#str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
775
776	if interpolation_points is not None:
777	str += ' Interpolation points (xi, eta):'\
778	' number of points == %d\n' %interpolation_points.shape[0]
779	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
780	max(interpolation_points[:,0]))
781	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
782	max(interpolation_points[:,1]))
783	str += ' Interpolated quantities (over all timesteps):\n'
784
785	for name in quantity_names:
786	q = precomputed_values[name][:].flat
787	str += ' %s at interpolation points in [%f, %f]\n'\
788	%(name, min(q), max(q))
789	str += '------------------------------------------------\n'
790
791	return str
792
793
794	def interpolate_sww(sww_file, time, interpolation_points,
795	quantity_names=None, verbose=False):
796	"""
797	obsolete.
798	use file_function in utils
799	"""
800	#open sww file
801	x, y, volumes, time, quantities = read_sww(sww_file)
802	print "x",x
803	print "y",y
804
805	print "time", time
806	print "quantities", quantities
807
808	#Add the x and y together
809	vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
810
811	#Will return the quantity values at the specified times and locations
812	interp = Interpolation_interface(time,
813	quantities,
814	quantity_names=quantity_names,
815	vertex_coordinates=vertex_coordinates,
816	triangles=volumes,
817	interpolation_points=interpolation_points,
818	verbose=verbose)
819
820
821	def read_sww(file_name):
822	"""
823	obsolete - Nothing should be calling this
824
825	Read in an sww file.
826
827	Input;
828	file_name - the sww file
829
830	Output;
831	x - Vector of x values
832	y - Vector of y values
833	z - Vector of bed elevation
834	volumes - Array. Each row has 3 values, representing
835	the vertices that define the volume
836	time - Vector of the times where there is stage information
837	stage - array with respect to time and vertices (x,y)
838	"""
839
840	msg = 'Function read_sww in interpolat.py is obsolete'
841	raise Exception, msg
842
843	#FIXME Have this reader as part of data_manager?
844
845	from Scientific.IO.NetCDF import NetCDFFile
846	import tempfile
847	import sys
848	import os
849
850	#Check contents
851	#Get NetCDF
852
853	# see if the file is there. Throw a QUIET IO error if it isn't
854	# I don't think I could implement the above
855
856	#throws prints to screen if file not present
857	junk = tempfile.mktemp(".txt")
858	fd = open(junk,'w')
859	stdout = sys.stdout
860	sys.stdout = fd
861	fid = NetCDFFile(file_name, 'r')
862	sys.stdout = stdout
863	fd.close()
864	#clean up
865	os.remove(junk)
866
867	# Get the variables
868	x = fid.variables['x'][:]
869	y = fid.variables['y'][:]
870	volumes = fid.variables['volumes'][:]
871	time = fid.variables['time'][:]
872
873	keys = fid.variables.keys()
874	keys.remove("x")
875	keys.remove("y")
876	keys.remove("volumes")
877	keys.remove("time")
878	#Turn NetCDF objects into Numeric arrays
879	quantities = {}
880	for name in keys:
881	quantities[name] = fid.variables[name][:]
882
883	fid.close()
884	return x, y, volumes, time, quantities
885
886
887	#-------------------------------------------------------------
888	if __name__ == "__main__":
889	names = ["x","y"]
890	someiterable = [[1,2],[3,4]]
891	csvwriter = writer(file("some.csv", "wb"))
892	csvwriter.writerow(names)
893	for row in someiterable:
894	csvwriter.writerow(row)

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

Download in other formats: