source: branches/numpy/anuga/fit_interpolate/interpolate.py @ 6689

Last change on this file since 6689 was 6689, checked in by rwilson, 15 years ago

Back-merge from Numeric trunk to numpy branch.

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