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

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

Refactored MeshQuad? into a self-contained class without global elements.

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