Context Navigation

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

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

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

Download in other formats: