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

Last change on this file since 6185 was 6185, checked in by ole, 15 years ago

Moved interpolate_polyline out of Interpolation class - it had nothing to do with it anyway.

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