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

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

Merged trunk into numpy, all tests and validations work.

File size: 47.7 KB
Line 
1"""Least squares interpolation.
2
3   Implements a least-squares interpolation.
4   Putting mesh data onto points.
5
6   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
7   Geoscience Australia, 2004.
8
9DESIGN ISSUES
10* what variables should be global?
11- if there are no global vars functions can be moved around alot easier
12
13* The public interface to Interpolate
14__init__
15interpolate
16interpolate_block
17
18"""
19
20import time
21import os
22import sys
23from warnings import warn
24from math import sqrt
25from csv import writer, DictWriter
26
27from anuga.caching.caching import cache
28from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
29from anuga.utilities.sparse import Sparse, Sparse_CSR
30from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError
31from anuga.coordinate_transforms.geo_reference import Geo_reference
32from anuga.utilities.numerical_tools import ensure_numeric, 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            if triangles is not None and vertex_coordinates is not None:
838                # Check that all interpolation points fall within
839                # mesh boundary as defined by triangles and vertex_coordinates.
840                from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
841                from anuga.utilities.polygon import outside_polygon
842
843                # Create temporary mesh object from mesh info passed
844                # into this function.
845                mesh = Mesh(vertex_coordinates, triangles)
846                mesh_boundary_polygon = mesh.get_boundary_polygon()
847
848                indices = outside_polygon(interpolation_points,
849                                          mesh_boundary_polygon)
850
851                # Record result
852                #self.mesh_boundary_polygon = mesh_boundary_polygon
853                self.indices_outside_mesh = indices
854
855                # Report
856                if len(indices) > 0:
857                    msg = 'Interpolation points in Interpolation function fall '
858                    msg += 'outside specified mesh. Offending points:\n'
859                    out_interp_pts = []
860                    for i in indices:
861                        msg += '%d: %s\n' % (i, interpolation_points[i])
862                        out_interp_pts.append(
863                                    ensure_numeric(interpolation_points[i]))
864
865                    if verbose is True:
866                        import sys
867                        if sys.platform == 'win32':
868                            # FIXME (Ole): Why only Windoze?
869                            from anuga.utilities.polygon import plot_polygons
870                            title = ('Interpolation points fall '
871                                     'outside specified mesh')
872                            plot_polygons([mesh_boundary_polygon,
873                                           interpolation_points,
874                                           out_interp_pts],
875                                          ['line', 'point', 'outside'],
876                                          figname='points_boundary_out',
877                                          label=title,
878                                          verbose=verbose)
879
880                    # Joaquim Luis suggested this as an Exception, so
881                    # that the user can now what the problem is rather than
882                    # looking for NaN's. However, NANs are handy as they can
883                    # be ignored leaving good points for continued processing.
884                    if verbose:
885                        print msg
886                    #raise Exception(msg)
887
888            elif triangles is None and vertex_coordinates is not None:    #jj
889                #Dealing with sts file
890                pass
891            else:
892                raise Exception('Sww file function requires both triangles and '
893                                'vertex_coordinates. sts file file function '
894                                'requires the latter.')
895
896            # Plot boundary and interpolation points
897            if verbose is True:
898                import sys
899                if sys.platform == 'win32':
900                    from anuga.utilities.polygon import plot_polygons
901                    title = ('Interpolation function: '
902                             'Polygon and interpolation points')
903                    plot_polygons([mesh_boundary_polygon,
904                                   interpolation_points],
905                                  ['line', 'point'],
906                                  figname='points_boundary',
907                                  label=title,
908                                  verbose=verbose)
909
910            m = len(self.interpolation_points)
911            p = len(self.time)
912
913            for name in quantity_names:
914                self.precomputed_values[name] = num.zeros((p, m), num.float)
915
916            if verbose is True:
917                print 'Build interpolator'
918
919
920            # Build interpolator
921            if triangles is not None and vertex_coordinates is not None:
922                if verbose:
923                    msg = 'Building interpolation matrix from source mesh '
924                    msg += '(%d vertices, %d triangles)' \
925                           % (vertex_coordinates.shape[0],
926                              triangles.shape[0])
927                    print msg
928
929                # This one is no longer needed for STS files
930                interpol = Interpolate(vertex_coordinates,
931                                       triangles,
932                                       verbose=verbose)
933
934            elif triangles is None and vertex_coordinates is not None:
935                if verbose:
936                    msg = 'Interpolation from STS file'
937                    print msg
938
939
940
941            if verbose:
942                print 'Interpolating (%d interpolation points, %d timesteps).' \
943                      %(self.interpolation_points.shape[0], self.time.shape[0]),
944
945                if time_thinning > 1:
946                    print 'Timesteps were thinned by a factor of %d' \
947                          % time_thinning
948                else:
949                    print
950
951            for i, t in enumerate(self.time):
952                # Interpolate quantities at this timestep
953                if verbose and i%((p+10)/10) == 0:
954                    print '  time step %d of %d' %(i, p)
955
956                for name in quantity_names:
957                    if len(quantities[name].shape) == 2:
958                        Q = quantities[name][i,:] # Quantities at timestep i
959                    else:
960                        Q = quantities[name][:]   # No time dependency
961
962                    if verbose and i%((p+10)/10) == 0:
963                        print '    quantity %s, size=%d' %(name, len(Q))
964
965                    # Interpolate
966                    if triangles is not None and vertex_coordinates is not None:
967                        result = interpol.interpolate(Q,
968                                                      point_coordinates=\
969                                                      self.interpolation_points,
970                                                      verbose=False) # No clutter
971                    elif triangles is None and vertex_coordinates is not None:
972                        result = interpolate_polyline(Q,
973                                                      vertex_coordinates,
974                                                      gauge_neighbour_id,
975                                                      interpolation_points=\
976                                                          self.interpolation_points)
977
978                    #assert len(result), len(interpolation_points)
979                    self.precomputed_values[name][i, :] = result
980
981            # Report
982            if verbose:
983                print self.statistics()
984        else:
985            # Store quantitites as is
986            for name in quantity_names:
987                self.precomputed_values[name] = quantities[name]
988
989    ##
990    # @brief Override object representation method.
991    def __repr__(self):
992        # return 'Interpolation function (spatio-temporal)'
993        return self.statistics()
994
995    ##
996    # @brief Override object() method.
997    # @param t Model time - must lie within existing timesteps.
998    # @param point_id Index of one of the preprocessed points.
999    # @param x ??
1000    # @param y ??
1001    # @return ??
1002    def __call__(self, t, point_id=None, x=None, y=None):
1003        """Evaluate f(t) or f(t, point_id)
1004
1005        Inputs:
1006          t:        time - Model time. Must lie within existing timesteps
1007          point_id: index of one of the preprocessed points.
1008
1009          If spatial info is present and all of point_id
1010          are None an exception is raised
1011
1012          If no spatial info is present, point_id arguments are ignored
1013          making f a function of time only.
1014
1015          FIXME: f(t, x, y) x, y could overrided location, point_id ignored
1016          FIXME: point_id could also be a slice
1017          FIXME: What if x and y are vectors?
1018          FIXME: What about f(x,y) without t?
1019        """
1020
1021        from math import pi, cos, sin, sqrt
1022        from anuga.abstract_2d_finite_volumes.util import mean
1023
1024        if self.spatial is True:
1025            if point_id is None:
1026                if x is None or y is None:
1027                    msg = 'Either point_id or x and y must be specified'
1028                    raise Exception(msg)
1029            else:
1030                if self.interpolation_points is None:
1031                    msg = 'Interpolation_function must be instantiated ' + \
1032                          'with a list of interpolation points before ' + \
1033                          'parameter point_id can be used'
1034                    raise Exception(msg)
1035
1036        msg = 'Time interval [%.16f:%.16f]' % (self.time[0], self.time[-1])
1037        msg += ' does not match model time: %.16f\n' % t
1038        if t < self.time[0]: raise Modeltime_too_early(msg)
1039        if t > self.time[-1]: raise Modeltime_too_late(msg)
1040
1041        oldindex = self.index #Time index
1042
1043        # Find current time slot
1044        while t > self.time[self.index]: self.index += 1
1045        while t < self.time[self.index]: self.index -= 1
1046
1047        if t == self.time[self.index]:
1048            # Protect against case where t == T[-1] (last time)
1049            #  - also works in general when t == T[i]
1050            ratio = 0
1051        else:
1052            # t is now between index and index+1
1053            ratio = ((t - self.time[self.index]) /
1054                         (self.time[self.index+1] - self.time[self.index]))
1055
1056        # Compute interpolated values
1057        q = num.zeros(len(self.quantity_names), num.float)
1058        for i, name in enumerate(self.quantity_names):
1059            Q = self.precomputed_values[name]
1060
1061            if self.spatial is False:
1062                # If there is no spatial info
1063                assert len(Q.shape) == 1
1064
1065                Q0 = Q[self.index]
1066                if ratio > 0: Q1 = Q[self.index+1]
1067            else:
1068                if x is not None and y is not None:
1069                    # Interpolate to x, y
1070                    raise 'x,y interpolation not yet implemented'
1071                else:
1072                    # Use precomputed point
1073                    Q0 = Q[self.index, point_id]
1074                    if ratio > 0:
1075                        Q1 = Q[self.index+1, point_id]
1076
1077            # Linear temporal interpolation
1078            if ratio > 0:
1079                if Q0 == NAN and Q1 == NAN:
1080                    q[i] = Q0
1081                else:
1082                    q[i] = Q0 + ratio*(Q1 - Q0)
1083            else:
1084                q[i] = Q0
1085
1086        # Return vector of interpolated values
1087        # FIXME:
1088        if self.spatial is True:
1089            return q
1090        else:
1091            # Replicate q according to x and y
1092            # This is e.g used for Wind_stress
1093            if x is None or y is None:
1094                return q
1095            else:
1096                try:
1097                    N = len(x)
1098                except:
1099                    return q
1100                else:
1101                    # x is a vector - Create one constant column for each value
1102                    N = len(x)
1103                    assert len(y) == N, 'x and y must have same length'
1104                    res = []
1105                    for col in q:
1106                        res.append(col*num.ones(N, num.float))
1107
1108                return res
1109
1110    ##
1111    # @brief Return model time as a vector of timesteps.
1112    def get_time(self):
1113        """Return model time as a vector of timesteps
1114        """
1115        return self.time
1116
1117    ##
1118    # @brief Output statistics about interpolation_function.
1119    # @return The statistics string.
1120    def statistics(self):
1121        """Output statistics about interpolation_function
1122        """
1123
1124        vertex_coordinates = self.vertex_coordinates
1125        interpolation_points = self.interpolation_points
1126        quantity_names = self.quantity_names
1127        #quantities = self.quantities
1128        precomputed_values = self.precomputed_values
1129
1130        x = vertex_coordinates[:,0]
1131        y = vertex_coordinates[:,1]
1132
1133        str =  '------------------------------------------------\n'
1134        str += 'Interpolation_function (spatio-temporal) statistics:\n'
1135        str += '  Extent:\n'
1136        str += '    x in [%f, %f], len(x) == %d\n'\
1137               %(min(x), max(x), len(x))
1138        str += '    y in [%f, %f], len(y) == %d\n'\
1139               %(min(y), max(y), len(y))
1140        str += '    t in [%f, %f], len(t) == %d\n'\
1141               %(min(self.time), max(self.time), len(self.time))
1142        str += '  Quantities:\n'
1143        for name in quantity_names:
1144            minq, maxq = self.quantities_range[name]
1145            str += '    %s in [%f, %f]\n' %(name, minq, maxq)
1146            #q = quantities[name][:].flatten()
1147            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
1148
1149        if interpolation_points is not None:
1150            str += '  Interpolation points (xi, eta):'\
1151                   ' number of points == %d\n' %interpolation_points.shape[0]
1152            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
1153                                            max(interpolation_points[:,0]))
1154            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
1155                                             max(interpolation_points[:,1]))
1156            str += '  Interpolated quantities (over all timesteps):\n'
1157
1158            for name in quantity_names:
1159                q = precomputed_values[name][:].flatten()
1160                str += '    %s at interpolation points in [%f, %f]\n'\
1161                       %(name, min(q), max(q))
1162        str += '------------------------------------------------\n'
1163
1164        return str
1165
1166
1167##
1168# @brief ??
1169# @param sww_file ??
1170# @param time ??
1171# @param interpolation_points ??
1172# @param quantity_names ??
1173# @param verbose ??
1174# @note Obsolete.  Use file_function() in utils.
1175def interpolate_sww(sww_file, time, interpolation_points,
1176                    quantity_names=None, verbose=False):
1177    """
1178    obsolete.
1179    use file_function in utils
1180    """
1181
1182    #open sww file
1183    x, y, volumes, time, quantities = read_sww(sww_file)
1184    print "x",x
1185    print "y",y
1186
1187    print "time", time
1188    print "quantities", quantities
1189
1190    #Add the x and y together
1191    vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),
1192                                         axis=1)
1193
1194    #Will return the quantity values at the specified times and locations
1195    interp = Interpolation_interface(time,
1196                                     quantities,
1197                                     quantity_names=quantity_names,
1198                                     vertex_coordinates=vertex_coordinates,
1199                                     triangles=volumes,
1200                                     interpolation_points=interpolation_points,
1201                                     verbose=verbose)
1202
1203
1204##
1205# @brief ??
1206# @param file_name Name of the .SWW file to read.
1207def read_sww(file_name):
1208    """
1209    obsolete - Nothing should be calling this
1210
1211    Read in an sww file.
1212
1213    Input;
1214    file_name - the sww file
1215
1216    Output;
1217    x - Vector of x values
1218    y - Vector of y values
1219    z - Vector of bed elevation
1220    volumes - Array.  Each row has 3 values, representing
1221    the vertices that define the volume
1222    time - Vector of the times where there is stage information
1223    stage - array with respect to time and vertices (x,y)
1224    """
1225
1226    msg = 'Function read_sww in interpolat.py is obsolete'
1227    raise Exception, msg
1228
1229    #FIXME Have this reader as part of data_manager?
1230
1231    from Scientific.IO.NetCDF import NetCDFFile
1232    import tempfile
1233    import sys
1234    import os
1235
1236    #Check contents
1237    #Get NetCDF
1238
1239    # see if the file is there.  Throw a QUIET IO error if it isn't
1240    # I don't think I could implement the above
1241
1242    #throws prints to screen if file not present
1243    junk = tempfile.mktemp(".txt")
1244    fd = open(junk,'w')
1245    stdout = sys.stdout
1246    sys.stdout = fd
1247    fid = NetCDFFile(file_name, netcdf_mode_r)
1248    sys.stdout = stdout
1249    fd.close()
1250    #clean up
1251    os.remove(junk)
1252
1253    # Get the variables
1254    x = fid.variables['x'][:]
1255    y = fid.variables['y'][:]
1256    volumes = fid.variables['volumes'][:]
1257    time = fid.variables['time'][:]
1258
1259    keys = fid.variables.keys()
1260    keys.remove("x")
1261    keys.remove("y")
1262    keys.remove("volumes")
1263    keys.remove("time")
1264     #Turn NetCDF objects into numeric arrays
1265    quantities = {}
1266    for name in keys:
1267        quantities[name] = fid.variables[name][:]
1268
1269    fid.close()
1270    return x, y, volumes, time, quantities
1271
1272
1273#-------------------------------------------------------------
1274if __name__ == "__main__":
1275    names = ["x","y"]
1276    someiterable = [[1,2],[3,4]]
1277    csvwriter = writer(file("some.csv", "wb"))
1278    csvwriter.writerow(names)
1279    for row in someiterable:
1280        csvwriter.writerow(row)
Note: See TracBrowser for help on using the repository browser.