Context Navigation

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

Last change on this file since 5657 was 5657, checked in by ole, 16 years ago
Implemented default_boundary option in File_boundary and Field_boundary as per ticket:293 and added a note in the user manual.
File size: 43.4 KB

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

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

Download in other formats: