Context Navigation

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

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

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

Download in other formats: