Context Navigation

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

Last change on this file since 4668 was 4663, checked in by duncan, 18 years ago
getting rid of xya more. Cleaning up.
File size: 32.2 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,
87	triangles,
88	mesh_origin,
89	verbose,
90	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	m = len(self.interpolation_points)
524	p = len(self.time)
525
526	for name in quantity_names:
527	self.precomputed_values[name] = zeros((p, m), Float)
528
529	# Build interpolator
530	if verbose:
531	msg = 'Building interpolation matrix from source mesh '
532	msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
533	triangles.shape[0])
534	print msg
535
536
537	interpol = Interpolate(vertex_coordinates,
538	triangles,
539	verbose=verbose)
540
541	if verbose:
542	print 'Interpolating (%d interpolation points, %d timesteps).'\
543	%(self.interpolation_points.shape[0], self.time.shape[0]),
544
545	if time_thinning > 1:
546	print 'Timesteps were thinned by a factor of %d' %time_thinning
547	else:
548	print
549
550	for i, t in enumerate(self.time):
551	# Interpolate quantities at this timestep
552	if verbose and i%((p+10)/10)==0:
553	print ' time step %d of %d' %(i, p)
554
555	for name in quantity_names:
556	if len(quantities[name].shape) == 2:
557	Q = quantities[name][i,:] # Quantities at timestep i
558	else:
559	Q = quantities[name][:] # No time dependency
560
561	# Interpolate
562	result = interpol.interpolate(Q,
563	point_coordinates=\
564	self.interpolation_points)
565
566	#assert len(result), len(interpolation_points)
567	self.precomputed_values[name][i, :] = result
568
569
570	# Report
571	if verbose:
572	print self.statistics()
573	#self.print_statistics()
574
575	else:
576	# Store quantitites as is
577	for name in quantity_names:
578	self.precomputed_values[name] = quantities[name]
579
580	def __repr__(self):
581	# return 'Interpolation function (spatio-temporal)'
582	return self.statistics()
583
584	def __call__(self, t, point_id=None, x=None, y=None):
585	"""Evaluate f(t) or f(t, point_id)
586
587	Inputs:
588	t: time - Model time. Must lie within existing timesteps
589	point_id: index of one of the preprocessed points.
590
591
592	If spatial info is present and all of point_id
593	are None an exception is raised
594
595	If no spatial info is present, point_id arguments are ignored
596	making f a function of time only.
597
598
599	FIXME: f(t, x, y) x, y could overrided location, point_id ignored
600	FIXME: point_id could also be a slice
601	FIXME: What if x and y are vectors?
602	FIXME: What about f(x,y) without t?
603	"""
604
605	from math import pi, cos, sin, sqrt
606	from Numeric import zeros, Float
607	from anuga.abstract_2d_finite_volumes.util import mean
608
609	if self.spatial is True:
610	if point_id is None:
611	if x is None or y is None:
612	msg = 'Either point_id or x and y must be specified'
613	raise Exception(msg)
614	else:
615	if self.interpolation_points is None:
616	msg = 'Interpolation_function must be instantiated ' +\
617	'with a list of interpolation points before parameter ' +\
618	'point_id can be used'
619	raise Exception(msg)
620
621	msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
622	msg += ' does not match model time: %.16f\n' %t
623	if t < self.time[0]: raise Exception(msg)
624	if t > self.time[-1]: raise Exception(msg)
625
626	oldindex = self.index #Time index
627
628	#Find current time slot
629	while t > self.time[self.index]: self.index += 1
630	while t < self.time[self.index]: self.index -= 1
631
632	if t == self.time[self.index]:
633	#Protect against case where t == T[-1] (last time)
634	# - also works in general when t == T[i]
635	ratio = 0
636	else:
637	#t is now between index and index+1
638	ratio = (t - self.time[self.index])/\
639	(self.time[self.index+1] - self.time[self.index])
640
641	#Compute interpolated values
642	q = zeros(len(self.quantity_names), Float)
643	#print "self.precomputed_values", self.precomputed_values
644	for i, name in enumerate(self.quantity_names):
645	Q = self.precomputed_values[name]
646
647	if self.spatial is False:
648	#If there is no spatial info
649	assert len(Q.shape) == 1
650
651	Q0 = Q[self.index]
652	if ratio > 0: Q1 = Q[self.index+1]
653
654	else:
655	if x is not None and y is not None:
656	#Interpolate to x, y
657
658	raise 'x,y interpolation not yet implemented'
659	else:
660	#Use precomputed point
661	Q0 = Q[self.index, point_id]
662	if ratio > 0:
663	Q1 = Q[self.index+1, point_id]
664
665	#Linear temporal interpolation
666	if ratio > 0:
667	if Q0 == NAN and Q1 == NAN:
668	q[i] = Q0
669	else:
670	q[i] = Q0 + ratio*(Q1 - Q0)
671	else:
672	q[i] = Q0
673
674
675	#Return vector of interpolated values
676	#if len(q) == 1:
677	# return q[0]
678	#else:
679	# return q
680
681
682	#Return vector of interpolated values
683	#FIXME:
684	if self.spatial is True:
685	return q
686	else:
687	#Replicate q according to x and y
688	#This is e.g used for Wind_stress
689	if x is None or y is None:
690	return q
691	else:
692	try:
693	N = len(x)
694	except:
695	return q
696	else:
697	from Numeric import ones, Float
698	#x is a vector - Create one constant column for each value
699	N = len(x)
700	assert len(y) == N, 'x and y must have same length'
701	res = []
702	for col in q:
703	res.append(col*ones(N, Float))
704
705	return res
706
707
708	def get_time(self):
709	"""Return model time as a vector of timesteps
710	"""
711	return self.time
712
713
714	def statistics(self):
715	"""Output statistics about interpolation_function
716	"""
717
718	vertex_coordinates = self.vertex_coordinates
719	interpolation_points = self.interpolation_points
720	quantity_names = self.quantity_names
721	#quantities = self.quantities
722	precomputed_values = self.precomputed_values
723
724	x = vertex_coordinates[:,0]
725	y = vertex_coordinates[:,1]
726
727	str = '------------------------------------------------\n'
728	str += 'Interpolation_function (spatio-temporal) statistics:\n'
729	str += ' Extent:\n'
730	str += ' x in [%f, %f], len(x) == %d\n'\
731	%(min(x), max(x), len(x))
732	str += ' y in [%f, %f], len(y) == %d\n'\
733	%(min(y), max(y), len(y))
734	str += ' t in [%f, %f], len(t) == %d\n'\
735	%(min(self.time), max(self.time), len(self.time))
736	str += ' Quantities:\n'
737	for name in quantity_names:
738	minq, maxq = self.quantities_range[name]
739	str += ' %s in [%f, %f]\n' %(name, minq, maxq)
740	#q = quantities[name][:].flat
741	#str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
742
743	if interpolation_points is not None:
744	str += ' Interpolation points (xi, eta):'\
745	' number of points == %d\n' %interpolation_points.shape[0]
746	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
747	max(interpolation_points[:,0]))
748	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
749	max(interpolation_points[:,1]))
750	str += ' Interpolated quantities (over all timesteps):\n'
751
752	for name in quantity_names:
753	q = precomputed_values[name][:].flat
754	str += ' %s at interpolation points in [%f, %f]\n'\
755	%(name, min(q), max(q))
756	str += '------------------------------------------------\n'
757
758	return str
759
760
761	def interpolate_sww(sww_file, time, interpolation_points,
762	quantity_names=None, verbose=False):
763	"""
764	obsolete.
765	use file_function in utils
766	"""
767	#open sww file
768	x, y, volumes, time, quantities = read_sww(sww_file)
769	print "x",x
770	print "y",y
771
772	print "time", time
773	print "quantities", quantities
774
775	#Add the x and y together
776	vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
777
778	#Will return the quantity values at the specified times and locations
779	interp = Interpolation_interface(time,
780	quantities,
781	quantity_names=quantity_names,
782	vertex_coordinates=vertex_coordinates,
783	triangles=volumes,
784	interpolation_points=interpolation_points,
785	verbose=verbose)
786
787
788	def read_sww(file_name):
789	"""
790	obsolete - Nothing should be calling this
791
792	Read in an sww file.
793
794	Input;
795	file_name - the sww file
796
797	Output;
798	x - Vector of x values
799	y - Vector of y values
800	z - Vector of bed elevation
801	volumes - Array. Each row has 3 values, representing
802	the vertices that define the volume
803	time - Vector of the times where there is stage information
804	stage - array with respect to time and vertices (x,y)
805	"""
806
807	msg = 'Function read_sww in interpolat.py is obsolete'
808	raise Exception, msg
809
810	#FIXME Have this reader as part of data_manager?
811
812	from Scientific.IO.NetCDF import NetCDFFile
813	import tempfile
814	import sys
815	import os
816
817	#Check contents
818	#Get NetCDF
819
820	# see if the file is there. Throw a QUIET IO error if it isn't
821	# I don't think I could implement the above
822
823	#throws prints to screen if file not present
824	junk = tempfile.mktemp(".txt")
825	fd = open(junk,'w')
826	stdout = sys.stdout
827	sys.stdout = fd
828	fid = NetCDFFile(file_name, 'r')
829	sys.stdout = stdout
830	fd.close()
831	#clean up
832	os.remove(junk)
833
834	# Get the variables
835	x = fid.variables['x'][:]
836	y = fid.variables['y'][:]
837	volumes = fid.variables['volumes'][:]
838	time = fid.variables['time'][:]
839
840	keys = fid.variables.keys()
841	keys.remove("x")
842	keys.remove("y")
843	keys.remove("volumes")
844	keys.remove("time")
845	#Turn NetCDF objects into Numeric arrays
846	quantities = {}
847	for name in keys:
848	quantities[name] = fid.variables[name][:]
849
850	fid.close()
851	return x, y, volumes, time, quantities
852
853
854	#-------------------------------------------------------------
855	if __name__ == "__main__":
856	names = ["x","y"]
857	someiterable = [[1,2],[3,4]]
858	csvwriter = writer(file("some.csv", "wb"))
859	csvwriter.writerow(names)
860	for row in someiterable:
861	csvwriter.writerow(row)

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

Download in other formats: