Context Navigation

source: branches/source_numpy_conversion/anuga/fit_interpolate/interpolate.py @ 6982

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

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

Download in other formats: