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

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

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