Context Navigation

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

Last change on this file since 4522 was 4522, checked in by duncan, 17 years ago
working on ticket#168
File size: 31.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.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	# for ticket 160
245	#boundary = self.mesh.get_boundary_polygon()
246	#geo = Geospatial_data(boundary)
247	#geo.export_points_file('serial-boundary.xya')
248	#geo.export_points_file('serial-boundary.txt')
249	if verbose: print 'Getting indices inside mesh boundary'
250	self.inside_poly_indices, self.outside_poly_indices = \
251	in_and_outside_polygon(point_coordinates,
252	self.mesh.get_boundary_polygon(),
253	closed = True, verbose = verbose)
254	#print "self.inside_poly_indices",self.inside_poly_indices
255	#print "self.outside_poly_indices",self.outside_poly_indices
256	#Build n x m interpolation matrix
257	if verbose and len(self.outside_poly_indices) >0:
258	print '\n WARNING: Points outside mesh boundary. \n'
259	# Since you can block, throw a warning, not an error.
260	if verbose and 0 == len(self.inside_poly_indices):
261	print '\n WARNING: No points within the mesh! \n'
262
263	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
264	n = point_coordinates.shape[0] # Nbr of data points
265
266	if verbose: print 'Number of datapoints: %d' %n
267	if verbose: print 'Number of basis functions: %d' %m
268
269	A = Sparse(n,m)
270
271	n = len(self.inside_poly_indices)
272	#Compute matrix elements for points inside the mesh
273	if verbose: print 'Building interpolation matrix fram %d points' %n
274	for k, i in enumerate(self.inside_poly_indices):
275	#For each data_coordinate point
276	if verbose and k%((n+10)/10)==0: print 'Doing %d of %d' %(k, n)
277	x = point_coordinates[i]
278	element_found, sigma0, sigma1, sigma2, k = \
279	search_tree_of_vertices(self.root, self.mesh, x)
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	self.precomputed_values[name][i, :] = result
551
552
553	# Report
554	if verbose:
555	print self.statistics()
556	#self.print_statistics()
557
558	else:
559	# Store quantitites as is
560	for name in quantity_names:
561	self.precomputed_values[name] = quantities[name]
562
563	def __repr__(self):
564	# return 'Interpolation function (spatio-temporal)'
565	return self.statistics()
566
567	def __call__(self, t, point_id=None, x=None, y=None):
568	"""Evaluate f(t) or f(t, point_id)
569
570	Inputs:
571	t: time - Model time. Must lie within existing timesteps
572	point_id: index of one of the preprocessed points.
573
574
575	If spatial info is present and all of point_id
576	are None an exception is raised
577
578	If no spatial info is present, point_id arguments are ignored
579	making f a function of time only.
580
581
582	FIXME: f(t, x, y) x, y could overrided location, point_id ignored
583	FIXME: point_id could also be a slice
584	FIXME: What if x and y are vectors?
585	FIXME: What about f(x,y) without t?
586	"""
587
588	from math import pi, cos, sin, sqrt
589	from Numeric import zeros, Float
590	from anuga.abstract_2d_finite_volumes.util import mean
591
592	if self.spatial is True:
593	if point_id is None:
594	if x is None or y is None:
595	msg = 'Either point_id or x and y must be specified'
596	raise Exception(msg)
597	else:
598	if self.interpolation_points is None:
599	msg = 'Interpolation_function must be instantiated ' +\
600	'with a list of interpolation points before parameter ' +\
601	'point_id can be used'
602	raise Exception(msg)
603
604	msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
605	msg += ' does not match model time: %.16f\n' %t
606	if t < self.time[0]: raise Exception(msg)
607	if t > self.time[-1]: raise Exception(msg)
608
609	oldindex = self.index #Time index
610
611	#Find current time slot
612	while t > self.time[self.index]: self.index += 1
613	while t < self.time[self.index]: self.index -= 1
614
615	if t == self.time[self.index]:
616	#Protect against case where t == T[-1] (last time)
617	# - also works in general when t == T[i]
618	ratio = 0
619	else:
620	#t is now between index and index+1
621	ratio = (t - self.time[self.index])/\
622	(self.time[self.index+1] - self.time[self.index])
623
624	#Compute interpolated values
625	q = zeros(len(self.quantity_names), Float)
626	#print "self.precomputed_values", self.precomputed_values
627	for i, name in enumerate(self.quantity_names):
628	Q = self.precomputed_values[name]
629
630	if self.spatial is False:
631	#If there is no spatial info
632	assert len(Q.shape) == 1
633
634	Q0 = Q[self.index]
635	if ratio > 0: Q1 = Q[self.index+1]
636
637	else:
638	if x is not None and y is not None:
639	#Interpolate to x, y
640
641	raise 'x,y interpolation not yet implemented'
642	else:
643	#Use precomputed point
644	Q0 = Q[self.index, point_id]
645	if ratio > 0:
646	Q1 = Q[self.index+1, point_id]
647
648	#Linear temporal interpolation
649	if ratio > 0:
650	if Q0 == NAN and Q1 == NAN:
651	q[i] = Q0
652	else:
653	q[i] = Q0 + ratio*(Q1 - Q0)
654	else:
655	q[i] = Q0
656
657
658	#Return vector of interpolated values
659	#if len(q) == 1:
660	# return q[0]
661	#else:
662	# return q
663
664
665	#Return vector of interpolated values
666	#FIXME:
667	if self.spatial is True:
668	return q
669	else:
670	#Replicate q according to x and y
671	#This is e.g used for Wind_stress
672	if x is None or y is None:
673	return q
674	else:
675	try:
676	N = len(x)
677	except:
678	return q
679	else:
680	from Numeric import ones, Float
681	#x is a vector - Create one constant column for each value
682	N = len(x)
683	assert len(y) == N, 'x and y must have same length'
684	res = []
685	for col in q:
686	res.append(col*ones(N, Float))
687
688	return res
689
690
691	def get_time(self):
692	"""Return model time as a vector of timesteps
693	"""
694	return self.time
695
696
697	def statistics(self):
698	"""Output statistics about interpolation_function
699	"""
700
701	vertex_coordinates = self.vertex_coordinates
702	interpolation_points = self.interpolation_points
703	quantity_names = self.quantity_names
704	#quantities = self.quantities
705	precomputed_values = self.precomputed_values
706
707	x = vertex_coordinates[:,0]
708	y = vertex_coordinates[:,1]
709
710	str = '------------------------------------------------\n'
711	str += 'Interpolation_function (spatio-temporal) statistics:\n'
712	str += ' Extent:\n'
713	str += ' x in [%f, %f], len(x) == %d\n'\
714	%(min(x), max(x), len(x))
715	str += ' y in [%f, %f], len(y) == %d\n'\
716	%(min(y), max(y), len(y))
717	str += ' t in [%f, %f], len(t) == %d\n'\
718	%(min(self.time), max(self.time), len(self.time))
719	str += ' Quantities:\n'
720	for name in quantity_names:
721	minq, maxq = self.quantities_range[name]
722	str += ' %s in [%f, %f]\n' %(name, minq, maxq)
723	#q = quantities[name][:].flat
724	#str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
725
726	if interpolation_points is not None:
727	str += ' Interpolation points (xi, eta):'\
728	' number of points == %d\n' %interpolation_points.shape[0]
729	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
730	max(interpolation_points[:,0]))
731	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
732	max(interpolation_points[:,1]))
733	str += ' Interpolated quantities (over all timesteps):\n'
734
735	for name in quantity_names:
736	q = precomputed_values[name][:].flat
737	str += ' %s at interpolation points in [%f, %f]\n'\
738	%(name, min(q), max(q))
739	str += '------------------------------------------------\n'
740
741	return str
742
743
744	def interpolate_sww(sww_file, time, interpolation_points,
745	quantity_names=None, verbose=False):
746	"""
747	obsolete.
748	use file_function in utils
749	"""
750	#open sww file
751	x, y, volumes, time, quantities = read_sww(sww_file)
752	print "x",x
753	print "y",y
754
755	print "time", time
756	print "quantities", quantities
757
758	#Add the x and y together
759	vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
760
761	#Will return the quantity values at the specified times and locations
762	interp = Interpolation_interface(time,
763	quantities,
764	quantity_names=quantity_names,
765	vertex_coordinates=vertex_coordinates,
766	triangles=volumes,
767	interpolation_points=interpolation_points,
768	verbose=verbose)
769
770
771	def read_sww(file_name):
772	"""
773	obsolete - Nothing should be calling this
774
775	Read in an sww file.
776
777	Input;
778	file_name - the sww file
779
780	Output;
781	x - Vector of x values
782	y - Vector of y values
783	z - Vector of bed elevation
784	volumes - Array. Each row has 3 values, representing
785	the vertices that define the volume
786	time - Vector of the times where there is stage information
787	stage - array with respect to time and vertices (x,y)
788	"""
789
790	#FIXME Have this reader as part of data_manager?
791
792	from Scientific.IO.NetCDF import NetCDFFile
793	import tempfile
794	import sys
795	import os
796
797	#Check contents
798	#Get NetCDF
799
800	# see if the file is there. Throw a QUIET IO error if it isn't
801	# I don't think I could implement the above
802
803	#throws prints to screen if file not present
804	junk = tempfile.mktemp(".txt")
805	fd = open(junk,'w')
806	stdout = sys.stdout
807	sys.stdout = fd
808	fid = NetCDFFile(file_name, 'r')
809	sys.stdout = stdout
810	fd.close()
811	#clean up
812	os.remove(junk)
813
814	# Get the variables
815	x = fid.variables['x'][:]
816	y = fid.variables['y'][:]
817	volumes = fid.variables['volumes'][:]
818	time = fid.variables['time'][:]
819
820	keys = fid.variables.keys()
821	keys.remove("x")
822	keys.remove("y")
823	keys.remove("volumes")
824	keys.remove("time")
825	#Turn NetCDF objects into Numeric arrays
826	quantities = {}
827	for name in keys:
828	quantities[name] = fid.variables[name][:]
829
830	fid.close()
831	return x, y, volumes, time, quantities
832
833
834	#-------------------------------------------------------------
835	if __name__ == "__main__":
836	names = ["x","y"]
837	someiterable = [[1,2],[3,4]]
838	csvwriter = writer(file("some.csv", "wb"))
839	csvwriter.writerow(names)
840	for row in someiterable:
841	csvwriter.writerow(row)

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

Download in other formats: