Context Navigation

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

Last change on this file since 4614 was 4614, checked in by duncan, 17 years ago
Centralising max vertices per cell.
File size: 31.7 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	# for ticket 160
244	#boundary = self.mesh.get_boundary_polygon()
245	#geo = Geospatial_data(boundary)
246	#geo.export_points_file('serial-boundary.xya')
247	#geo.export_points_file('serial-boundary.txt')
248	if verbose: print 'Getting indices inside mesh boundary'
249	self.inside_poly_indices, self.outside_poly_indices = \
250	in_and_outside_polygon(point_coordinates,
251	self.mesh.get_boundary_polygon(),
252	closed = True, verbose = verbose)
253	#print "self.inside_poly_indices",self.inside_poly_indices
254	#print "self.outside_poly_indices",self.outside_poly_indices
255	#Build n x m interpolation matrix
256	if verbose and len(self.outside_poly_indices) > 0:
257	print '\n WARNING: Points outside mesh boundary. \n'
258	# Since you can block, throw a warning, not an error.
259	if verbose and 0 == len(self.inside_poly_indices):
260	print '\n WARNING: No points within the mesh! \n'
261
262	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
263	n = point_coordinates.shape[0] # Nbr of data points
264
265	if verbose: print 'Number of datapoints: %d' %n
266	if verbose: print 'Number of basis functions: %d' %m
267
268	A = Sparse(n,m)
269
270	n = len(self.inside_poly_indices)
271	#Compute matrix elements for points inside the mesh
272	if verbose: print 'Building interpolation matrix from %d points' %n
273	for d, i in enumerate(self.inside_poly_indices):
274	# For each data_coordinate point
275	if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
276	x = point_coordinates[i]
277	element_found, sigma0, sigma1, sigma2, k = \
278	search_tree_of_vertices(self.root, self.mesh, x)
279
280	# Update interpolation matrix A if necessary
281	if element_found is True:
282	# Assign values to matrix A
283
284	j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0
285	j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
286	j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
287
288	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
289	js = [j0,j1,j2]
290
291	for j in js:
292	A[i,j] = sigmas[j]
293	else:
294	msg = 'Could not find triangle for point', x
295	raise Exception(msg)
296	return A
297
298	def interpolate_sww2csv(sww_file,
299	points,
300	depth_file,
301	velocity_x_file,
302	velocity_y_file,
303	#quantities = ['depth', 'velocity'],
304	verbose=True,
305	use_cache = True):
306	"""
307	Interpolate the quantities at a given set of locations, given
308	an sww file.
309	The results are written to a csv file.
310
311	In the future let points be a csv or xya file.
312	And the user choose the quantities.
313
314	This is currently quite specific.
315	If it need to be more general, chagne things.
316
317	This is really returning speed, not velocity.
318	"""
319	quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
320	#print "points",points
321	points = ensure_absolute(points)
322	point_count = len(points)
323	callable_sww = file_function(sww_file,
324	quantities=quantities,
325	interpolation_points=points,
326	verbose=verbose,
327	use_cache=use_cache)
328
329	depth_writer = writer(file(depth_file, "wb"))
330	velocity_x_writer = writer(file(velocity_x_file, "wb"))
331	velocity_y_writer = writer(file(velocity_y_file, "wb"))
332	# Write heading
333	heading = [str(x[0])+ ':' + str(x[1]) for x in points]
334	heading.insert(0, "time")
335	depth_writer.writerow(heading)
336	velocity_x_writer.writerow(heading)
337	velocity_y_writer.writerow(heading)
338
339	for time in callable_sww.get_time():
340	depths = [time]
341	velocity_xs = [time]
342	velocity_ys = [time]
343	for point_i, point in enumerate(points):
344	quantities = callable_sww(time,point_i)
345	#print "quantities", quantities
346
347	w = quantities[0]
348	z = quantities[1]
349	momentum_x = quantities[2]
350	momentum_y = quantities[3]
351	depth = w - z
352
353	if w == NAN or z == NAN or momentum_x == NAN:
354	velocity_x = NAN
355	else:
356	if depth > 1.e-30: # use epsilon
357	velocity_x = momentum_x / depth #Absolute velocity
358	else:
359	velocity_x = 0
360	if w == NAN or z == NAN or momentum_y == NAN:
361	velocity_y = NAN
362	else:
363	if depth > 1.e-30: # use epsilon
364	velocity_y = momentum_y / depth #Absolute velocity
365	else:
366	velocity_y = 0
367	depths.append(depth)
368	velocity_xs.append(velocity_x)
369	velocity_ys.append(velocity_y)
370	depth_writer.writerow(depths)
371	velocity_x_writer.writerow(velocity_xs)
372	velocity_y_writer.writerow(velocity_ys)
373
374
375	class Interpolation_function:
376	"""Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
377	which is interpolated from time series defined at vertices of
378	triangular mesh (such as those stored in sww files)
379
380	Let m be the number of vertices, n the number of triangles
381	and p the number of timesteps.
382
383	Mandatory input
384	time: px1 array of monotonously increasing times (Float)
385	quantities: Dictionary of arrays or 1 array (Float)
386	The arrays must either have dimensions pxm or mx1.
387	The resulting function will be time dependent in
388	the former case while it will be constan with
389	respect to time in the latter case.
390
391	Optional input:
392	quantity_names: List of keys into the quantities dictionary
393	vertex_coordinates: mx2 array of coordinates (Float)
394	triangles: nx3 array of indices into vertex_coordinates (Int)
395	interpolation_points: Nx2 array of coordinates to be interpolated to
396	verbose: Level of reporting
397
398
399	The quantities returned by the callable object are specified by
400	the list quantities which must contain the names of the
401	quantities to be returned and also reflect the order, e.g. for
402	the shallow water wave equation, on would have
403	quantities = ['stage', 'xmomentum', 'ymomentum']
404
405	The parameter interpolation_points decides at which points interpolated
406	quantities are to be computed whenever object is called.
407	If None, return average value
408
409	FIXME (Ole): Need to allow vertex coordinates and interpolation points to be
410	geospatial data objects
411
412	Time assumed to be relative to starttime
413	"""
414
415
416	def __init__(self,
417	time,
418	quantities,
419	quantity_names=None,
420	vertex_coordinates=None,
421	triangles=None,
422	interpolation_points=None,
423	time_thinning=1,
424	verbose=False):
425	"""Initialise object and build spatial interpolation if required
426
427	Time_thinning_number controls how many timesteps to use. Only timesteps with
428	index%time_thinning_number == 0 will used, or in other words a value of 3, say,
429	will cause the algorithm to use every third time step.
430	"""
431
432	from Numeric import array, zeros, Float, alltrue, concatenate,\
433	reshape, ArrayType
434
435
436	from anuga.config import time_format
437	import types
438
439
440	# Check temporal info
441	time = ensure_numeric(time)
442	msg = 'Time must be a monotonuosly '
443	msg += 'increasing sequence %s' %time
444	assert alltrue(time[1:] - time[:-1] >= 0 ), msg
445
446
447	# Check if quantities is a single array only
448	if type(quantities) != types.DictType:
449	quantities = ensure_numeric(quantities)
450	quantity_names = ['Attribute']
451
452	#Make it a dictionary
453	quantities = {quantity_names[0]: quantities}
454
455
456	# Use keys if no names are specified
457	if quantity_names is None:
458	quantity_names = quantities.keys()
459
460
461	# Check spatial info
462	if vertex_coordinates is None:
463	self.spatial = False
464	else:
465	vertex_coordinates = ensure_absolute(vertex_coordinates)
466
467	assert triangles is not None, 'Triangles array must be specified'
468	triangles = ensure_numeric(triangles)
469	self.spatial = True
470
471	# Thin timesteps if needed
472	# Note array() is used to make the thinned arrays contiguous in memory
473	self.time = array(time[::time_thinning])
474	for name in quantity_names:
475	if len(quantities[name].shape) == 2:
476	quantities[name] = array(quantities[name][::time_thinning,:])
477
478	# Save for use with statistics
479	self.quantities_range = {}
480	for name in quantity_names:
481	q = quantities[name][:].flat
482	self.quantities_range[name] = [min(q), max(q)]
483
484	self.quantity_names = quantity_names
485	self.vertex_coordinates = vertex_coordinates
486	self.interpolation_points = interpolation_points
487
488
489	self.index = 0 # Initial time index
490	self.precomputed_values = {}
491
492
493	# Precomputed spatial interpolation if requested
494	if interpolation_points is not None:
495	if self.spatial is False:
496	raise 'Triangles and vertex_coordinates must be specified'
497
498	try:
499	self.interpolation_points = ensure_numeric(interpolation_points)
500	except:
501	msg = 'Interpolation points must be an N x 2 Numeric array '+\
502	'or a list of points\n'
503	msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
504	'...')
505	raise msg
506
507
508	m = len(self.interpolation_points)
509	p = len(self.time)
510
511	for name in quantity_names:
512	self.precomputed_values[name] = zeros((p, m), Float)
513
514	# Build interpolator
515	if verbose:
516	msg = 'Building interpolation matrix from source mesh '
517	msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
518	triangles.shape[0])
519	print msg
520
521
522	interpol = Interpolate(vertex_coordinates,
523	triangles,
524	verbose=verbose)
525
526	if verbose:
527	print 'Interpolating (%d interpolation points, %d timesteps).'\
528	%(self.interpolation_points.shape[0], self.time.shape[0]),
529
530	if time_thinning > 1:
531	print 'Timesteps were thinned by a factor of %d' %time_thinning
532	else:
533	print
534
535	for i, t in enumerate(self.time):
536	# Interpolate quantities at this timestep
537	if verbose and i%((p+10)/10)==0:
538	print ' time step %d of %d' %(i, p)
539
540	for name in quantity_names:
541	if len(quantities[name].shape) == 2:
542	Q = quantities[name][i,:] # Quantities at timestep i
543	else:
544	Q = quantities[name][:] # No time dependency
545
546	# Interpolate
547	result = interpol.interpolate(Q,
548	point_coordinates=\
549	self.interpolation_points)
550
551	#assert len(result), len(interpolation_points)
552	self.precomputed_values[name][i, :] = result
553
554
555	# Report
556	if verbose:
557	print self.statistics()
558	#self.print_statistics()
559
560	else:
561	# Store quantitites as is
562	for name in quantity_names:
563	self.precomputed_values[name] = quantities[name]
564
565	def __repr__(self):
566	# return 'Interpolation function (spatio-temporal)'
567	return self.statistics()
568
569	def __call__(self, t, point_id=None, x=None, y=None):
570	"""Evaluate f(t) or f(t, point_id)
571
572	Inputs:
573	t: time - Model time. Must lie within existing timesteps
574	point_id: index of one of the preprocessed points.
575
576
577	If spatial info is present and all of point_id
578	are None an exception is raised
579
580	If no spatial info is present, point_id arguments are ignored
581	making f a function of time only.
582
583
584	FIXME: f(t, x, y) x, y could overrided location, point_id ignored
585	FIXME: point_id could also be a slice
586	FIXME: What if x and y are vectors?
587	FIXME: What about f(x,y) without t?
588	"""
589
590	from math import pi, cos, sin, sqrt
591	from Numeric import zeros, Float
592	from anuga.abstract_2d_finite_volumes.util import mean
593
594	if self.spatial is True:
595	if point_id is None:
596	if x is None or y is None:
597	msg = 'Either point_id or x and y must be specified'
598	raise Exception(msg)
599	else:
600	if self.interpolation_points is None:
601	msg = 'Interpolation_function must be instantiated ' +\
602	'with a list of interpolation points before parameter ' +\
603	'point_id can be used'
604	raise Exception(msg)
605
606	msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
607	msg += ' does not match model time: %.16f\n' %t
608	if t < self.time[0]: raise Exception(msg)
609	if t > self.time[-1]: raise Exception(msg)
610
611	oldindex = self.index #Time index
612
613	#Find current time slot
614	while t > self.time[self.index]: self.index += 1
615	while t < self.time[self.index]: self.index -= 1
616
617	if t == self.time[self.index]:
618	#Protect against case where t == T[-1] (last time)
619	# - also works in general when t == T[i]
620	ratio = 0
621	else:
622	#t is now between index and index+1
623	ratio = (t - self.time[self.index])/\
624	(self.time[self.index+1] - self.time[self.index])
625
626	#Compute interpolated values
627	q = zeros(len(self.quantity_names), Float)
628	#print "self.precomputed_values", self.precomputed_values
629	for i, name in enumerate(self.quantity_names):
630	Q = self.precomputed_values[name]
631
632	if self.spatial is False:
633	#If there is no spatial info
634	assert len(Q.shape) == 1
635
636	Q0 = Q[self.index]
637	if ratio > 0: Q1 = Q[self.index+1]
638
639	else:
640	if x is not None and y is not None:
641	#Interpolate to x, y
642
643	raise 'x,y interpolation not yet implemented'
644	else:
645	#Use precomputed point
646	Q0 = Q[self.index, point_id]
647	if ratio > 0:
648	Q1 = Q[self.index+1, point_id]
649
650	#Linear temporal interpolation
651	if ratio > 0:
652	if Q0 == NAN and Q1 == NAN:
653	q[i] = Q0
654	else:
655	q[i] = Q0 + ratio*(Q1 - Q0)
656	else:
657	q[i] = Q0
658
659
660	#Return vector of interpolated values
661	#if len(q) == 1:
662	# return q[0]
663	#else:
664	# return q
665
666
667	#Return vector of interpolated values
668	#FIXME:
669	if self.spatial is True:
670	return q
671	else:
672	#Replicate q according to x and y
673	#This is e.g used for Wind_stress
674	if x is None or y is None:
675	return q
676	else:
677	try:
678	N = len(x)
679	except:
680	return q
681	else:
682	from Numeric import ones, Float
683	#x is a vector - Create one constant column for each value
684	N = len(x)
685	assert len(y) == N, 'x and y must have same length'
686	res = []
687	for col in q:
688	res.append(col*ones(N, Float))
689
690	return res
691
692
693	def get_time(self):
694	"""Return model time as a vector of timesteps
695	"""
696	return self.time
697
698
699	def statistics(self):
700	"""Output statistics about interpolation_function
701	"""
702
703	vertex_coordinates = self.vertex_coordinates
704	interpolation_points = self.interpolation_points
705	quantity_names = self.quantity_names
706	#quantities = self.quantities
707	precomputed_values = self.precomputed_values
708
709	x = vertex_coordinates[:,0]
710	y = vertex_coordinates[:,1]
711
712	str = '------------------------------------------------\n'
713	str += 'Interpolation_function (spatio-temporal) statistics:\n'
714	str += ' Extent:\n'
715	str += ' x in [%f, %f], len(x) == %d\n'\
716	%(min(x), max(x), len(x))
717	str += ' y in [%f, %f], len(y) == %d\n'\
718	%(min(y), max(y), len(y))
719	str += ' t in [%f, %f], len(t) == %d\n'\
720	%(min(self.time), max(self.time), len(self.time))
721	str += ' Quantities:\n'
722	for name in quantity_names:
723	minq, maxq = self.quantities_range[name]
724	str += ' %s in [%f, %f]\n' %(name, minq, maxq)
725	#q = quantities[name][:].flat
726	#str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
727
728	if interpolation_points is not None:
729	str += ' Interpolation points (xi, eta):'\
730	' number of points == %d\n' %interpolation_points.shape[0]
731	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
732	max(interpolation_points[:,0]))
733	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
734	max(interpolation_points[:,1]))
735	str += ' Interpolated quantities (over all timesteps):\n'
736
737	for name in quantity_names:
738	q = precomputed_values[name][:].flat
739	str += ' %s at interpolation points in [%f, %f]\n'\
740	%(name, min(q), max(q))
741	str += '------------------------------------------------\n'
742
743	return str
744
745
746	def interpolate_sww(sww_file, time, interpolation_points,
747	quantity_names=None, verbose=False):
748	"""
749	obsolete.
750	use file_function in utils
751	"""
752	#open sww file
753	x, y, volumes, time, quantities = read_sww(sww_file)
754	print "x",x
755	print "y",y
756
757	print "time", time
758	print "quantities", quantities
759
760	#Add the x and y together
761	vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
762
763	#Will return the quantity values at the specified times and locations
764	interp = Interpolation_interface(time,
765	quantities,
766	quantity_names=quantity_names,
767	vertex_coordinates=vertex_coordinates,
768	triangles=volumes,
769	interpolation_points=interpolation_points,
770	verbose=verbose)
771
772
773	def read_sww(file_name):
774	"""
775	obsolete - Nothing should be calling this
776
777	Read in an sww file.
778
779	Input;
780	file_name - the sww file
781
782	Output;
783	x - Vector of x values
784	y - Vector of y values
785	z - Vector of bed elevation
786	volumes - Array. Each row has 3 values, representing
787	the vertices that define the volume
788	time - Vector of the times where there is stage information
789	stage - array with respect to time and vertices (x,y)
790	"""
791
792	msg = 'Function read_sww in interpolat.py is obsolete'
793	raise Exception, msg
794
795	#FIXME Have this reader as part of data_manager?
796
797	from Scientific.IO.NetCDF import NetCDFFile
798	import tempfile
799	import sys
800	import os
801
802	#Check contents
803	#Get NetCDF
804
805	# see if the file is there. Throw a QUIET IO error if it isn't
806	# I don't think I could implement the above
807
808	#throws prints to screen if file not present
809	junk = tempfile.mktemp(".txt")
810	fd = open(junk,'w')
811	stdout = sys.stdout
812	sys.stdout = fd
813	fid = NetCDFFile(file_name, 'r')
814	sys.stdout = stdout
815	fd.close()
816	#clean up
817	os.remove(junk)
818
819	# Get the variables
820	x = fid.variables['x'][:]
821	y = fid.variables['y'][:]
822	volumes = fid.variables['volumes'][:]
823	time = fid.variables['time'][:]
824
825	keys = fid.variables.keys()
826	keys.remove("x")
827	keys.remove("y")
828	keys.remove("volumes")
829	keys.remove("time")
830	#Turn NetCDF objects into Numeric arrays
831	quantities = {}
832	for name in keys:
833	quantities[name] = fid.variables[name][:]
834
835	fid.close()
836	return x, y, volumes, time, quantities
837
838
839	#-------------------------------------------------------------
840	if __name__ == "__main__":
841	names = ["x","y"]
842	someiterable = [[1,2],[3,4]]
843	csvwriter = writer(file("some.csv", "wb"))
844	csvwriter.writerow(names)
845	for row in someiterable:
846	csvwriter.writerow(row)

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

Download in other formats: