Context Navigation

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

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

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

Download in other formats: