Context Navigation

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

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

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

Download in other formats: