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

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

Ticket 344 completed and unit tests pass.
Stripped out a heap of tabs and replaced them with 4 spaces.
Removed some redundant imports and deprecated functions.

File size: 50.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, 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 interpolated f
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        # Quick test against boundary, but will not deal with holes in the mesh
467        inside_boundary_indices, outside_poly_indices = \
468            in_and_outside_polygon(point_coordinates,
469                                   self.mesh.get_boundary_polygon(),
470                                   closed=True, verbose=verbose)
471
472        # Build n x m interpolation matrix
473        if verbose and len(outside_poly_indices) > 0:
474            log.critical('WARNING: Points outside mesh boundary.')
475
476        # Since you can block, throw a warning, not an error.
477        if verbose and 0 == len(inside_boundary_indices):
478            log.critical('WARNING: No points within the mesh!')
479
480        m = self.mesh.number_of_nodes  # Nbr of basis functions (1/vertex)
481        n = point_coordinates.shape[0] # Nbr of data points
482
483        if verbose: log.critical('Number of datapoints: %d' % n)
484        if verbose: log.critical('Number of basis functions: %d' % m)
485
486        A = Sparse(n,m)
487
488        n = len(inside_boundary_indices)
489
490        centroids = []
491       
492        inside_poly_indices = []
493       
494        # Compute matrix elements for points inside the mesh
495        if verbose: log.critical('Building interpolation matrix from %d points'
496                                 % n)
497
498        for d, i in enumerate(inside_boundary_indices):
499            # For each data_coordinate point
500            if verbose and d%((n+10)/10)==0: log.critical('Doing %d of %d'
501                                                          %(d, n))
502
503            x = point_coordinates[i]
504            element_found, sigma0, sigma1, sigma2, k = \
505                           search_tree_of_vertices(self.root, x)
506                       
507        # Update interpolation matrix A if necessary
508            if element_found is True:
509
510                if verbose:
511                    print 'Point is within mesh:', d, i           
512           
513                inside_poly_indices.append(i)
514               
515                # Assign values to matrix A
516                j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0
517                j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
518                j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
519                js = [j0, j1, j2]
520               
521                if output_centroids is False:
522                    # Weight each vertex according to its distance from x
523                    sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
524                    for j in js:
525                        A[i, j] = sigmas[j]
526                else:
527                    # If centroids are needed, weight all 3 vertices equally
528                    for j in js:
529                        A[i, j] = 1.0/3.0
530                    centroids.append(self.mesh.centroid_coordinates[k])                       
531            else:
532                if verbose:
533                    log.critical('Mesh has a hole - moving this point to outside list')
534
535                # This is a numpy arrays, so we need to do a slow transfer
536                outside_poly_indices = num.append(outside_poly_indices, [i], axis=0)
537
538        return A, inside_poly_indices, outside_poly_indices, centroids
539
540
541
542
543
544
545##
546# @brief ??
547# @param vertices ??
548# @param vertex_attributes ??
549# @param triangles ??
550# @param points ??
551# @param max_points_per_cell ??
552# @param start_blocking_len ??
553# @param mesh_origin ??
554def benchmark_interpolate(vertices,
555                          vertex_attributes,
556                          triangles, points,
557                          max_points_per_cell=None,
558                          start_blocking_len=500000,
559                          mesh_origin=None):
560    """
561    points: Interpolate mesh data to these positions.
562            List of coordinate pairs [x, y] of
563            data points or an nx2 numeric array or a Geospatial_data object
564
565    No test for this yet.
566    Note, this has no time the input data has no time dimension.  Which is
567    different from most of the data we interpolate, eg sww info.
568
569    Output:
570        Interpolated values at inputted points.
571    """
572
573    interp = Interpolate(vertices,
574                         triangles,
575                         max_vertices_per_cell=max_points_per_cell,
576                         mesh_origin=mesh_origin)
577
578    calc = interp.interpolate(vertex_attributes,
579                              points,
580                              start_blocking_len=start_blocking_len)
581
582
583##
584# @brief Interpolate quantities at given locations (from .SWW file).
585# @param sww_file Input .SWW file.
586# @param points A list of the 'gauges' x,y location.
587# @param depth_file The name of the output depth file.
588# @param velocity_x_file Name of the output x velocity  file.
589# @param velocity_y_file Name of the output y velocity  file.
590# @param stage_file Name of the output stage file.
591# @param froude_file
592# @param time_thinning Time thinning step to use.
593# @param verbose True if this function is to be verbose.
594# @param use_cache True if we are caching.
595def interpolate_sww2csv(sww_file,
596                        points,
597                        depth_file,
598                        velocity_x_file,
599                        velocity_y_file,
600                        stage_file=None,
601                        froude_file=None,
602                        time_thinning=1,
603                        verbose=True,
604                        use_cache = True):
605    """
606    Interpolate the quantities at a given set of locations, given
607    an sww file.
608    The results are written to csv files.
609
610    sww_file is the input sww file.
611    points is a list of the 'gauges' x,y location.
612    depth_file is the name of the output depth file
613    velocity_x_file is the name of the output x velocity file.
614    velocity_y_file is the name of the output y velocity file.
615    stage_file is the name of the output stage file.
616
617    In the csv files columns represents the gauges and each row is a
618    time slice.
619
620    Time_thinning_number controls how many timesteps to use. Only
621    timesteps with index%time_thinning_number == 0 will used, or
622    in other words a value of 3, say, will cause the algorithm to
623    use every third time step.
624
625    In the future let points be a points file.
626    And let the user choose the quantities.
627
628    This is currently quite specific.
629    If it is need to be more general, change things.
630    """
631
632    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
633    points = ensure_absolute(points)
634    point_count = len(points)
635    callable_sww = file_function(sww_file,
636                                 quantities=quantities,
637                                 interpolation_points=points,
638                                 verbose=verbose,
639                                 time_thinning=time_thinning,
640                                 use_cache=use_cache)
641
642    depth_writer = writer(file(depth_file, "wb"))
643    velocity_x_writer = writer(file(velocity_x_file, "wb"))
644    velocity_y_writer = writer(file(velocity_y_file, "wb"))
645    if stage_file is not None:
646        stage_writer = writer(file(stage_file, "wb"))
647    if froude_file is not None:
648        froude_writer = writer(file(froude_file, "wb"))
649
650    # Write heading
651    heading = [str(x[0])+ ':' + str(x[1]) for x in points]
652    heading.insert(0, "time")
653    depth_writer.writerow(heading)
654    velocity_x_writer.writerow(heading)
655    velocity_y_writer.writerow(heading)
656    if stage_file is not None:
657        stage_writer.writerow(heading)
658    if froude_file is not None:
659        froude_writer.writerow(heading)
660
661    for time in callable_sww.get_time():
662        depths = [time]
663        velocity_xs = [time]
664        velocity_ys = [time]
665        if stage_file is not None:
666            stages = [time]
667        if froude_file is not None:
668            froudes = [time]
669        for point_i, point in enumerate(points):
670            quantities = callable_sww(time,point_i)
671
672            w = quantities[0]
673            z = quantities[1]
674            momentum_x = quantities[2]
675            momentum_y = quantities[3]
676            depth = w - z
677
678            if w == NAN or z == NAN or momentum_x == NAN:
679                velocity_x = NAN
680            else:
681                if depth > 1.e-30: # use epsilon
682                    velocity_x = momentum_x / depth  #Absolute velocity
683                else:
684                    velocity_x = 0
685
686            if w == NAN or z == NAN or momentum_y == NAN:
687                velocity_y = NAN
688            else:
689                if depth > 1.e-30: # use epsilon
690                    velocity_y = momentum_y / depth  #Absolute velocity
691                else:
692                    velocity_y = 0
693
694            if depth < 1.e-30: # use epsilon
695                froude = NAN
696            else:
697                froude = sqrt(velocity_x*velocity_x + velocity_y*velocity_y)/ \
698                         sqrt(depth * 9.8066) # gravity m/s/s
699
700            depths.append(depth)
701            velocity_xs.append(velocity_x)
702            velocity_ys.append(velocity_y)
703
704            if stage_file is not None:
705                stages.append(w)
706            if froude_file is not None:
707                froudes.append(froude)
708
709        depth_writer.writerow(depths)
710        velocity_x_writer.writerow(velocity_xs)
711        velocity_y_writer.writerow(velocity_ys)
712
713        if stage_file is not None:
714            stage_writer.writerow(stages)
715        if froude_file is not None:
716            froude_writer.writerow(froudes)
717
718
719##
720# @brief
721class Interpolation_function:
722    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
723    which is interpolated from time series defined at vertices of
724    triangular mesh (such as those stored in sww files)
725
726    Let m be the number of vertices, n the number of triangles
727    and p the number of timesteps.
728    Also, let N be the number of interpolation points.
729
730    Mandatory input
731        time:                 px1 array of monotonously increasing times (float)
732        quantities:           Dictionary of arrays or 1 array (float)
733                              The arrays must either have dimensions pxm or mx1.
734                              The resulting function will be time dependent in
735                              the former case while it will be constant with
736                              respect to time in the latter case.
737
738    Optional input:
739        quantity_names:       List of keys into the quantities dictionary for
740                              imposing a particular order on the output vector.
741        vertex_coordinates:   mx2 array of coordinates (float)
742        triangles:            nx3 array of indices into vertex_coordinates (int)
743        interpolation_points: Nx2 array of coordinates to be interpolated to
744        verbose:              Level of reporting
745
746    The quantities returned by the callable object are specified by
747    the list quantities which must contain the names of the
748    quantities to be returned and also reflect the order, e.g. for
749    the shallow water wave equation, on would have
750    quantities = ['stage', 'xmomentum', 'ymomentum']
751
752    The parameter interpolation_points decides at which points interpolated
753    quantities are to be computed whenever object is called.
754    If None, return average value
755
756    FIXME (Ole): Need to allow vertex coordinates and interpolation points to
757                 be geospatial data objects
758
759    (FIXME (Ole): This comment should be removed)
760    Time assumed to be relative to starttime
761    All coordinates assume origin of (0,0) - e.g. georeferencing must be
762    taken care of outside this function
763    """
764
765    ##
766    # @brief ??
767    # @param time ??
768    # @param quantities ??
769    # @param quantity_names   ??
770    # @param vertex_coordinates ??
771    # @param triangles ??
772    # @param interpolation_points ??
773    # @param time_thinning ??
774    # @param verbose ??
775    # @param gauge_neighbour_id ??
776    def __init__(self,
777                 time,
778                 quantities,
779                 quantity_names=None,
780                 vertex_coordinates=None,
781                 triangles=None,
782                 interpolation_points=None,
783                 time_thinning=1,
784                 verbose=False,
785                 gauge_neighbour_id=None,
786                 output_centroids=False):
787        """Initialise object and build spatial interpolation if required
788
789        Time_thinning_number controls how many timesteps to use. Only timesteps
790        with index%time_thinning_number == 0 will used, or in other words a
791        value of 3, say, will cause the algorithm to use every third time step.
792        """
793
794        from anuga.config import time_format
795        import types
796
797        if verbose is True:
798            log.critical('Interpolation_function: input checks')
799
800        # Check temporal info
801        time = ensure_numeric(time)
802
803        if not num.alltrue(time[1:] - time[:-1] >= 0):
804            # This message is time consuming to form due to the conversion of
805            msg = 'Time must be a monotonuosly increasing sequence %s' % time
806            raise Exception, msg
807
808        # Check if quantities is a single array only
809        if type(quantities) != types.DictType:
810            quantities = ensure_numeric(quantities)
811            quantity_names = ['Attribute']
812
813            # Make it a dictionary
814            quantities = {quantity_names[0]: quantities}
815
816        # Use keys if no names are specified
817        if quantity_names is None:
818            quantity_names = quantities.keys()
819
820        # Check spatial info
821        if vertex_coordinates is None:
822            self.spatial = False
823        else:
824            # FIXME (Ole): Try ensure_numeric here -
825            #              this function knows nothing about georefering.
826            vertex_coordinates = ensure_absolute(vertex_coordinates)
827
828            if triangles is not None:
829                triangles = ensure_numeric(triangles)
830            self.spatial = True
831
832        if verbose is True:
833            log.critical('Interpolation_function: thinning by %d'
834                         % time_thinning)
835
836
837        # Thin timesteps if needed
838        # Note array() is used to make the thinned arrays contiguous in memory
839        self.time = num.array(time[::time_thinning])
840        for name in quantity_names:
841            if len(quantities[name].shape) == 2:
842                quantities[name] = num.array(quantities[name][::time_thinning,:])
843
844        if verbose is True:
845            log.critical('Interpolation_function: precomputing')
846
847        # Save for use with statistics
848        self.quantities_range = {}
849        for name in quantity_names:
850            q = quantities[name][:].flatten()
851            self.quantities_range[name] = [min(q), max(q)]
852
853        self.quantity_names = quantity_names
854        self.vertex_coordinates = vertex_coordinates
855        self.interpolation_points = interpolation_points
856
857        self.index = 0    # Initial time index
858        self.precomputed_values = {}
859        self.centroids = []
860
861        # Precomputed spatial interpolation if requested
862        if interpolation_points is not None:
863            #no longer true. sts files have spatial = True but
864            #if self.spatial is False:
865            #    raise 'Triangles and vertex_coordinates must be specified'
866            #
867            try:
868                self.interpolation_points = \
869                    interpolation_points = ensure_numeric(interpolation_points)
870            except:
871                msg = 'Interpolation points must be an N x 2 numeric array ' \
872                      'or a list of points\n'
873                msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...')
874                raise msg
875
876            # Ensure 'mesh_boundary_polygon' is defined
877            mesh_boundary_polygon = None
878           
879            if triangles is not None and vertex_coordinates is not None:
880                # Check that all interpolation points fall within
881                # mesh boundary as defined by triangles and vertex_coordinates.
882                from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
883                from anuga.utilities.polygon import outside_polygon
884
885                # Create temporary mesh object from mesh info passed
886                # into this function.
887                mesh = Mesh(vertex_coordinates, triangles)
888                mesh_boundary_polygon = mesh.get_boundary_polygon()
889
890                indices = outside_polygon(interpolation_points,
891                                          mesh_boundary_polygon)
892
893                # Record result
894                #self.mesh_boundary_polygon = mesh_boundary_polygon
895                self.indices_outside_mesh = indices
896
897                # Report
898                if len(indices) > 0:
899                    msg = 'Interpolation points in Interpolation function fall '
900                    msg += 'outside specified mesh. Offending points:\n'
901                    out_interp_pts = []
902                    for i in indices:
903                        msg += '%d: %s\n' % (i, interpolation_points[i])
904                        out_interp_pts.append(
905                                    ensure_numeric(interpolation_points[i]))
906
907                    if verbose is True:
908                        import sys
909                        if sys.platform == 'win32':
910                            # FIXME (Ole): Why only Windoze?
911                            from anuga.utilities.polygon import plot_polygons
912                            title = ('Interpolation points fall '
913                                     'outside specified mesh')
914                            plot_polygons([mesh_boundary_polygon,
915                                           interpolation_points,
916                                           out_interp_pts],
917                                          ['line', 'point', 'outside'],
918                                          figname='points_boundary_out',
919                                          label=title,
920                                          verbose=verbose)
921
922                    # Joaquim Luis suggested this as an Exception, so
923                    # that the user can now what the problem is rather than
924                    # looking for NaN's. However, NANs are handy as they can
925                    # be ignored leaving good points for continued processing.
926                    if verbose:
927                        log.critical(msg)
928                    #raise Exception(msg)
929
930            elif triangles is None and vertex_coordinates is not None:    #jj
931                #Dealing with sts file
932                pass
933            else:
934                raise Exception('Sww file function requires both triangles and '
935                                'vertex_coordinates. sts file file function '
936                                'requires the latter.')
937
938            # Plot boundary and interpolation points,
939            # but only if if 'mesh_boundary_polygon' has data.
940            if verbose is True and mesh_boundary_polygon is not None:
941                import sys
942                if sys.platform == 'win32':
943                    from anuga.utilities.polygon import plot_polygons
944                    title = ('Interpolation function: '
945                             'Polygon and interpolation points')
946                    plot_polygons([mesh_boundary_polygon,
947                                   interpolation_points],
948                                  ['line', 'point'],
949                                  figname='points_boundary',
950                                  label=title,
951                                  verbose=verbose)
952
953            m = len(self.interpolation_points)
954            p = len(self.time)
955
956            for name in quantity_names:
957                self.precomputed_values[name] = num.zeros((p, m), num.float)
958
959            if verbose is True:
960                log.critical('Build interpolator')
961
962
963            # Build interpolator
964            if triangles is not None and vertex_coordinates is not None:
965                if verbose:
966                    msg = 'Building interpolation matrix from source mesh '
967                    msg += '(%d vertices, %d triangles)' \
968                           % (vertex_coordinates.shape[0],
969                              triangles.shape[0])
970                    log.critical(msg)
971
972                # This one is no longer needed for STS files
973                interpol = Interpolate(vertex_coordinates,
974                                       triangles,
975                                       verbose=verbose)
976
977            elif triangles is None and vertex_coordinates is not None:
978                if verbose:
979                    log.critical('Interpolation from STS file')
980
981
982
983            if verbose:
984                log.critical('Interpolating (%d interpolation points, %d timesteps).'
985                             % (self.interpolation_points.shape[0], self.time.shape[0]))
986
987                if time_thinning > 1:
988                    log.critical('Timesteps were thinned by a factor of %d'
989                                 % time_thinning)
990                else:
991                    log.critical()
992
993            for i, t in enumerate(self.time):
994                # Interpolate quantities at this timestep
995                if verbose and i%((p+10)/10) == 0:
996                    log.critical('  time step %d of %d' % (i, p))
997
998                for name in quantity_names:
999                    if len(quantities[name].shape) == 2:
1000                        Q = quantities[name][i,:] # Quantities at timestep i
1001                    else:
1002                        Q = quantities[name][:]   # No time dependency
1003
1004                    if verbose and i%((p+10)/10) == 0:
1005                        log.critical('    quantity %s, size=%d' % (name, len(Q)))
1006
1007                    # Interpolate
1008                    if triangles is not None and vertex_coordinates is not None:
1009                        result = interpol.interpolate(Q,
1010                                                      point_coordinates=\
1011                                                      self.interpolation_points,
1012                                                      verbose=False,
1013                                                      output_centroids=output_centroids)
1014                        self.centroids = interpol.centroids                                                         
1015                    elif triangles is None and vertex_coordinates is not None:
1016                        result = interpolate_polyline(Q,
1017                                                      vertex_coordinates,
1018                                                      gauge_neighbour_id,
1019                                                      interpolation_points=\
1020                                                          self.interpolation_points)
1021
1022                    #assert len(result), len(interpolation_points)
1023                    self.precomputed_values[name][i, :] = result                                   
1024                   
1025            # Report
1026            if verbose:
1027                log.critical(self.statistics())           
1028        else:
1029            # Store quantitites as is
1030            for name in quantity_names:
1031                self.precomputed_values[name] = quantities[name]
1032
1033    ##
1034    # @brief Override object representation method.
1035    def __repr__(self):
1036        # return 'Interpolation function (spatio-temporal)'
1037        return self.statistics()
1038
1039    ##
1040    # @brief Evaluate interpolation function
1041    # @param t Model time - must lie within existing timesteps.
1042    # @param point_id Index of one of the preprocessed points.
1043    # @param x ??
1044    # @param y ??
1045    # @return ??
1046    def __call__(self, t, point_id=None, x=None, y=None):
1047        """Evaluate f(t) or f(t, point_id)
1048
1049        Inputs:
1050          t:        time - Model time. Must lie within existing timesteps
1051          point_id: index of one of the preprocessed points.
1052
1053          If spatial info is present and all of point_id
1054          are None an exception is raised
1055
1056          If no spatial info is present, point_id arguments are ignored
1057          making f a function of time only.
1058
1059          FIXME: f(t, x, y) x, y could overrided location, point_id ignored
1060          FIXME: point_id could also be a slice
1061          FIXME: What if x and y are vectors?
1062          FIXME: What about f(x,y) without t?
1063        """
1064
1065        from math import pi, cos, sin, sqrt
1066
1067        if self.spatial is True:
1068            if point_id is None:
1069                if x is None or y is None:
1070                    msg = 'Either point_id or x and y must be specified'
1071                    raise Exception(msg)
1072            else:
1073                if self.interpolation_points is None:
1074                    msg = 'Interpolation_function must be instantiated ' + \
1075                          'with a list of interpolation points before ' + \
1076                          'parameter point_id can be used'
1077                    raise Exception(msg)
1078
1079        msg = 'Time interval [%.16f:%.16f]' % (self.time[0], self.time[-1])
1080        msg += ' does not match model time: %.16f\n' % t
1081        if t < self.time[0]: raise Modeltime_too_early(msg)
1082        if t > self.time[-1]: raise Modeltime_too_late(msg)
1083
1084        oldindex = self.index #Time index
1085
1086        # Find current time slot
1087        while t > self.time[self.index]: self.index += 1
1088        while t < self.time[self.index]: self.index -= 1
1089
1090        if t == self.time[self.index]:
1091            # Protect against case where t == T[-1] (last time)
1092            #  - also works in general when t == T[i]
1093            ratio = 0
1094        else:
1095            # t is now between index and index+1
1096            ratio = ((t - self.time[self.index]) /
1097                         (self.time[self.index+1] - self.time[self.index]))
1098
1099        # Compute interpolated values
1100        q = num.zeros(len(self.quantity_names), num.float)
1101        for i, name in enumerate(self.quantity_names):
1102            Q = self.precomputed_values[name]
1103
1104            if self.spatial is False:
1105                # If there is no spatial info
1106                assert len(Q.shape) == 1
1107
1108                Q0 = Q[self.index]
1109                if ratio > 0: Q1 = Q[self.index+1]
1110            else:
1111                if x is not None and y is not None:
1112                    # Interpolate to x, y
1113                    raise 'x,y interpolation not yet implemented'
1114                else:
1115                    # Use precomputed point
1116                    Q0 = Q[self.index, point_id]
1117                    if ratio > 0:
1118                        Q1 = Q[self.index+1, point_id]
1119
1120            # Linear temporal interpolation
1121            if ratio > 0:
1122                if Q0 == NAN and Q1 == NAN:
1123                    q[i] = Q0
1124                else:
1125                    q[i] = Q0 + ratio*(Q1 - Q0)
1126            else:
1127                q[i] = Q0
1128
1129        # Return vector of interpolated values
1130        # FIXME:
1131        if self.spatial is True:
1132            return q
1133        else:
1134            # Replicate q according to x and y
1135            # This is e.g used for Wind_stress
1136            if x is None or y is None:
1137                return q
1138            else:
1139                try:
1140                    N = len(x)
1141                except:
1142                    return q
1143                else:
1144                    # x is a vector - Create one constant column for each value
1145                    N = len(x)
1146                    assert len(y) == N, 'x and y must have same length'
1147                    res = []
1148                    for col in q:
1149                        res.append(col*num.ones(N, num.float))
1150
1151                return res
1152
1153    ##
1154    # @brief Return model time as a vector of timesteps.
1155    def get_time(self):
1156        """Return model time as a vector of timesteps
1157        """
1158        return self.time
1159
1160    ##
1161    # @brief Output statistics about interpolation_function.
1162    # @return The statistics string.
1163    def statistics(self):
1164        """Output statistics about interpolation_function
1165        """
1166
1167        vertex_coordinates = self.vertex_coordinates
1168        interpolation_points = self.interpolation_points
1169        quantity_names = self.quantity_names
1170        #quantities = self.quantities
1171        precomputed_values = self.precomputed_values
1172
1173        x = vertex_coordinates[:,0]
1174        y = vertex_coordinates[:,1]
1175
1176        str =  '------------------------------------------------\n'
1177        str += 'Interpolation_function (spatio-temporal) statistics:\n'
1178        str += '  Extent:\n'
1179        str += '    x in [%f, %f], len(x) == %d\n'\
1180               %(min(x), max(x), len(x))
1181        str += '    y in [%f, %f], len(y) == %d\n'\
1182               %(min(y), max(y), len(y))
1183        str += '    t in [%f, %f], len(t) == %d\n'\
1184               %(min(self.time), max(self.time), len(self.time))
1185        str += '  Quantities:\n'
1186        for name in quantity_names:
1187            minq, maxq = self.quantities_range[name]
1188            str += '    %s in [%f, %f]\n' %(name, minq, maxq)
1189            #q = quantities[name][:].flatten()
1190            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
1191
1192        if interpolation_points is not None:
1193            str += '  Interpolation points (xi, eta):'\
1194                   ' number of points == %d\n' %interpolation_points.shape[0]
1195            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
1196                                            max(interpolation_points[:,0]))
1197            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
1198                                             max(interpolation_points[:,1]))
1199            str += '  Interpolated quantities (over all timesteps):\n'
1200
1201            for name in quantity_names:
1202                q = precomputed_values[name][:].flatten()
1203                str += '    %s at interpolation points in [%f, %f]\n'\
1204                       %(name, min(q), max(q))
1205        str += '------------------------------------------------\n'
1206
1207        return str
1208
1209
1210##
1211# @brief ??
1212# @param sww_file ??
1213# @param time ??
1214# @param interpolation_points ??
1215# @param quantity_names ??
1216# @param verbose ??
1217# @note Obsolete.  Use file_function() in utils.
1218def interpolate_sww(sww_file, time, interpolation_points,
1219                    quantity_names=None, verbose=False):
1220    """
1221    obsolete.
1222    use file_function in utils
1223    """
1224
1225    #open sww file
1226    x, y, volumes, time, quantities = read_sww(sww_file)
1227    log.critical("x=%s" % str(x))
1228    log.critical("y=%s" % str(y))
1229
1230    log.critical("time=%s" % str(time))
1231    log.critical("quantities=%s" % str(quantities))
1232
1233    #Add the x and y together
1234    vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),
1235                                         axis=1)
1236
1237    #Will return the quantity values at the specified times and locations
1238    interp = Interpolation_interface(time,
1239                                     quantities,
1240                                     quantity_names=quantity_names,
1241                                     vertex_coordinates=vertex_coordinates,
1242                                     triangles=volumes,
1243                                     interpolation_points=interpolation_points,
1244                                     verbose=verbose)
1245
1246
1247##
1248# @brief ??
1249# @param file_name Name of the .SWW file to read.
1250def read_sww(file_name):
1251    """
1252    obsolete - Nothing should be calling this
1253
1254    Read in an sww file.
1255
1256    Input;
1257    file_name - the sww file
1258
1259    Output;
1260    x - Vector of x values
1261    y - Vector of y values
1262    z - Vector of bed elevation
1263    volumes - Array.  Each row has 3 values, representing
1264    the vertices that define the volume
1265    time - Vector of the times where there is stage information
1266    stage - array with respect to time and vertices (x,y)
1267    """
1268
1269    msg = 'Function read_sww in interpolat.py is obsolete'
1270    raise Exception, msg
1271
1272    #FIXME Have this reader as part of data_manager?
1273
1274    from Scientific.IO.NetCDF import NetCDFFile
1275    import tempfile
1276    import sys
1277    import os
1278
1279    #Check contents
1280    #Get NetCDF
1281
1282    # see if the file is there.  Throw a QUIET IO error if it isn't
1283    # I don't think I could implement the above
1284
1285    #throws prints to screen if file not present
1286    junk = tempfile.mktemp(".txt")
1287    fd = open(junk,'w')
1288    stdout = sys.stdout
1289    sys.stdout = fd
1290    fid = NetCDFFile(file_name, netcdf_mode_r)
1291    sys.stdout = stdout
1292    fd.close()
1293    #clean up
1294    os.remove(junk)
1295
1296    # Get the variables
1297    x = fid.variables['x'][:]
1298    y = fid.variables['y'][:]
1299    volumes = fid.variables['volumes'][:]
1300    time = fid.variables['time'][:]
1301
1302    keys = fid.variables.keys()
1303    keys.remove("x")
1304    keys.remove("y")
1305    keys.remove("volumes")
1306    keys.remove("time")
1307     #Turn NetCDF objects into numeric arrays
1308    quantities = {}
1309    for name in keys:
1310        quantities[name] = fid.variables[name][:]
1311
1312    fid.close()
1313    return x, y, volumes, time, quantities
1314
1315
1316#-------------------------------------------------------------
1317if __name__ == "__main__":
1318    names = ["x","y"]
1319    someiterable = [[1,2],[3,4]]
1320    csvwriter = writer(file("some.csv", "wb"))
1321    csvwriter.writerow(names)
1322    for row in someiterable:
1323        csvwriter.writerow(row)
Note: See TracBrowser for help on using the repository browser.