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

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

Ticket 113 - added an output_centroids parameter to allow centroid data to be written to gauges, rather than the data at the sampled point.

Added 2 new unit tests for this functionality.

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