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

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

Added documentation. Fixes in sww2timeseries.

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