source: trunk/anuga_core/source/anuga/fit_interpolate/interpolate.py @ 8050

Last change on this file since 8050 was 7854, checked in by hudson, 14 years ago

Fixed Windows Interpolate nonexistent verbose flag.

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