Context Navigation

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

Last change on this file since 5361 was 5361, checked in by jakeman, 17 years ago
Adding ability to read sts file and create file function
File size: 41.3 KB

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

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

Download in other formats: