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

Last change on this file since 5196 was 5196, checked in by duncan, 16 years ago

Additional functionality. Output stage .csv as well.

File size: 35.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
22from warnings import warn
23from math import sqrt
24from csv import writer, DictWriter
25
26from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
27     ArrayType, allclose, take, NewAxis, arange
28
29from anuga.caching.caching import cache
30from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
31from anuga.utilities.sparse import Sparse, Sparse_CSR
32from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError
33from anuga.coordinate_transforms.geo_reference import Geo_reference
34from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
35from anuga.utilities.polygon import in_and_outside_polygon
36from anuga.geospatial_data.geospatial_data import Geospatial_data
37from anuga.geospatial_data.geospatial_data import ensure_absolute
38from anuga.fit_interpolate.search_functions import search_tree_of_vertices
39from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
40from anuga.abstract_2d_finite_volumes.util import file_function
41
42
43class Interpolate (FitInterpolate):
44       
45    def __init__(self,
46                 vertex_coordinates,
47                 triangles,
48                 mesh_origin=None,
49                 verbose=False,
50                 max_vertices_per_cell=None):
51
52
53        """ Build interpolation matrix mapping from
54        function values at vertices to function values at data points
55
56        Inputs:
57
58          vertex_coordinates: List of coordinate pairs [xi, eta] of
59              points constituting a mesh (or an m x 2 Numeric array or
60              a geospatial object)
61              Points may appear multiple times
62              (e.g. if vertices have discontinuities)
63
64          triangles: List of 3-tuples (or a Numeric array) of
65              integers representing indices of all vertices in the mesh.
66
67          mesh_origin: A geo_reference object or 3-tuples consisting of
68              UTM zone, easting and northing.
69              If specified vertex coordinates are assumed to be
70              relative to their respective origins.
71
72          max_vertices_per_cell: Number of vertices in a quad tree cell
73          at which the cell is split into 4.
74
75          Note: Don't supply a vertex coords as a geospatial object and
76              a mesh origin, since geospatial has its own mesh origin.
77        """
78
79        # FIXME (Ole): Need an input check
80       
81        # Initialise variabels
82        self._A_can_be_reused = False
83        self._point_coordinates = None
84       
85        FitInterpolate.__init__(self,
86                                vertex_coordinates=vertex_coordinates,
87                                triangles=triangles,
88                                mesh_origin=mesh_origin,
89                                verbose=verbose,
90                                max_vertices_per_cell=max_vertices_per_cell)
91
92    # FIXME: What is a good start_blocking_len value?
93    def interpolate(self,
94                    f,
95                    point_coordinates=None,
96                    start_blocking_len=500000,
97                    verbose=False):
98        """Interpolate mesh data f to determine values, z, at points.
99
100        f is the data on the mesh vertices.
101
102        The mesh values representing a smooth surface are
103        assumed to be specified in f.
104
105        Inputs:
106          f: Vector or array of data at the mesh vertices.
107              If f is an array, interpolation will be done for each column as
108              per underlying matrix-matrix multiplication
109          point_coordinates: Interpolate mesh data to these positions.
110              List of coordinate pairs [x, y] of
111              data points or an nx2 Numeric array or a Geospatial_data object
112             
113              If point_coordinates is absent, the points inputted last time
114              this method was called are used, if possible.
115          start_blocking_len: If the # of points is more or greater than this,
116              start blocking
117
118        Output:
119          Interpolated values at inputted points (z).
120        """
121
122       
123        # FIXME (Ole): Need an input check that dimensions are compatible
124
125        # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
126        # method is called even if interpolation points are unchanged.
127       
128        #print "point_coordinates interpolate.interpolate", point_coordinates
129        if isinstance(point_coordinates, Geospatial_data):
130            point_coordinates = point_coordinates.get_data_points( \
131                absolute = True)
132
133        # Can I interpolate, based on previous point_coordinates?
134        if point_coordinates is None:
135            if self._A_can_be_reused is True and \
136                   len(self._point_coordinates) < start_blocking_len:
137                z = self._get_point_data_z(f,
138                                           verbose=verbose)
139            elif self._point_coordinates is not None:
140                #     if verbose, give warning
141                if verbose:
142                    print 'WARNING: Recalculating A matrix, due to blocking.'
143                point_coordinates = self._point_coordinates
144            else:
145                #There are no good point_coordinates. import sys; sys.exit()
146                msg = 'ERROR (interpolate.py): No point_coordinates inputted'
147                raise Exception(msg)
148       
149        if point_coordinates is not None:   
150            self._point_coordinates = point_coordinates
151            if len(point_coordinates) < start_blocking_len or \
152                   start_blocking_len == 0:
153                self._A_can_be_reused = True
154                z = self.interpolate_block(f, point_coordinates,
155                                           verbose=verbose)
156            else:
157                #print 'BLOCKING'
158                #Handle blocking
159                self._A_can_be_reused = False
160                start = 0
161                # creating a dummy array to concatenate to.
162               
163                f = ensure_numeric(f, Float)
164                #print "f.shape",f.shape
165                if len(f.shape) > 1:
166                    z = zeros((0,f.shape[1]))
167                else:
168                    z = zeros((0,))
169                   
170                for end in range(start_blocking_len,
171                                 len(point_coordinates),
172                                 start_blocking_len):
173                   
174                    t = self.interpolate_block(f, point_coordinates[start:end],
175                                               verbose=verbose)
176                    #print "t", t
177                    #print "z", z
178                    z = concatenate((z,t))
179                    start = end
180                   
181                end = len(point_coordinates)
182                t = self.interpolate_block(f, point_coordinates[start:end],
183                                           verbose=verbose)
184                z = concatenate((z,t))
185        return z
186
187    def interpolate_block(self, f, point_coordinates=None, verbose=False):
188        """
189        Call this if you want to control the blocking or make sure blocking
190        doesn't occur.
191
192        Return the point data, z.
193       
194        See interpolate for doc info.
195        """
196        if isinstance(point_coordinates,Geospatial_data):
197            point_coordinates = point_coordinates.get_data_points( \
198                absolute = True)
199        if point_coordinates is not None:
200            self._A = self._build_interpolation_matrix_A(point_coordinates,
201                                                         verbose=verbose)
202        return self._get_point_data_z(f)
203
204    def _get_point_data_z(self, f, verbose=False):
205        """
206        Return the point data, z.
207       
208        Precondition,
209        The _A matrix has been created
210        """
211        z = self._A * f
212        # Taking into account points outside the mesh.
213        #print "self.outside_poly_indices", self.outside_poly_indices
214        #print "self.inside_poly_indices", self.inside_poly_indices
215        #print "z", z
216        for i in self.outside_poly_indices: 
217            z[i] = NAN
218        return z
219
220    def _build_interpolation_matrix_A(self,
221                                      point_coordinates,
222                                      verbose=False):
223        """Build n x m interpolation matrix, where
224        n is the number of data points and
225        m is the number of basis functions phi_k (one per vertex)
226
227        This algorithm uses a quad tree data structure for fast binning
228        of data points
229        origin is a 3-tuple consisting of UTM zone, easting and northing.
230        If specified coordinates are assumed to be relative to this origin.
231
232        This one will override any data_origin that may be specified in
233        instance interpolation
234
235        Preconditions
236        Point_coordindates and mesh vertices have the same origin.
237        """
238
239        #print 'Building interpolation matrix'
240       
241        #Convert point_coordinates to Numeric arrays, in case it was a list.
242        point_coordinates = ensure_numeric(point_coordinates, Float)
243       
244        if verbose: print 'Getting indices inside mesh boundary'
245        self.inside_poly_indices, self.outside_poly_indices  = \
246                     in_and_outside_polygon(point_coordinates,
247                                            self.mesh.get_boundary_polygon(),
248                                            closed = True, verbose = verbose)
249       
250        #Build n x m interpolation matrix
251        if verbose and len(self.outside_poly_indices) > 0:
252            print '\n WARNING: Points outside mesh boundary. \n'
253        # Since you can block, throw a warning, not an error.
254        if verbose and 0 == len(self.inside_poly_indices):
255            print '\n WARNING: No points within the mesh! \n'
256           
257        m = self.mesh.number_of_nodes  # Nbr of basis functions (1/vertex)
258        n = point_coordinates.shape[0] # Nbr of data points
259
260        if verbose: print 'Number of datapoints: %d' %n
261        if verbose: print 'Number of basis functions: %d' %m
262
263        A = Sparse(n,m)
264
265        n = len(self.inside_poly_indices)
266        #Compute matrix elements for points inside the mesh
267        if verbose: print 'Building interpolation matrix from %d points' %n
268        for d, i in enumerate(self.inside_poly_indices):
269            # For each data_coordinate point
270            if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
271            x = point_coordinates[i]
272            element_found, sigma0, sigma1, sigma2, k = \
273                           search_tree_of_vertices(self.root, self.mesh, x)
274           
275            # Update interpolation matrix A if necessary
276            if element_found is True:
277                # Assign values to matrix A
278
279                j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0
280                j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
281                j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
282
283                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
284                js     = [j0,j1,j2]
285
286                for j in js:
287                    A[i,j] = sigmas[j]
288            else:
289                msg = 'Could not find triangle for point', x
290                raise Exception(msg)
291        return A
292
293def benchmark_interpolate(vertices,
294                          vertex_attributes,
295                          triangles, points,
296                          max_points_per_cell=None,
297                          start_blocking_len=500000,
298                          mesh_origin=None):
299    """
300    points: Interpolate mesh data to these positions.
301              List of coordinate pairs [x, y] of
302              data points or an nx2 Numeric array or a Geospatial_data object
303             
304    No test for this yet.
305    Note, this has no time the input data has no time dimension.  Which is
306    different from most of the data we interpolate, eg sww info.
307     
308        Output:
309          Interpolated values at inputted points.
310    """
311    interp = Interpolate(vertices,
312                         triangles, 
313                         max_vertices_per_cell=max_points_per_cell,
314                         mesh_origin=mesh_origin)
315           
316    calc = interp.interpolate(vertex_attributes
317                              ,points
318                              ,start_blocking_len=start_blocking_len)
319    #print "calc", calc
320   
321def interpolate_sww2csv(sww_file,
322                        points,
323                        depth_file,
324                        velocity_x_file,
325                        velocity_y_file,
326                        stage_file=None,
327                        #quantities = ['depth', 'velocity'],
328                        verbose=True,
329                        use_cache = True):
330    """
331    Interpolate the quantities at a given set of locations, given
332    an sww file.
333    The results are written to a csv file.
334
335    In the future let points be a points file.
336    And the user choose the quantities.
337
338    This is currently quite specific.
339    If it need to be more general, chagne things.
340
341    This is really returning speed, not velocity.
342    """
343    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
344    #print "points",points
345    points = ensure_absolute(points)
346    point_count = len(points)
347    callable_sww = file_function(sww_file,
348                                 quantities=quantities,
349                                 interpolation_points=points,
350                                 verbose=verbose,
351                                 use_cache=use_cache)
352   
353    depth_writer = writer(file(depth_file, "wb"))
354    velocity_x_writer = writer(file(velocity_x_file, "wb"))
355    velocity_y_writer = writer(file(velocity_y_file, "wb"))
356    if stage_file is not None:
357        stage_writer = writer(file(stage_file, "wb"))
358    # Write heading
359    heading = [str(x[0])+ ':' + str(x[1]) for x in points]
360    heading.insert(0, "time")
361    depth_writer.writerow(heading)
362    velocity_x_writer.writerow(heading)
363    velocity_y_writer.writerow(heading)
364    if stage_file is not None:
365        stage_writer.writerow(heading)   
366   
367    for time in callable_sww.get_time():
368        depths = [time]
369        velocity_xs = [time]
370        velocity_ys = [time]
371        if stage_file is not None: 
372            stages = [time] 
373        for point_i, point in enumerate(points):
374            quantities = callable_sww(time,point_i)
375            #print "quantities", quantities
376           
377            w = quantities[0]
378            z = quantities[1]
379            momentum_x = quantities[2]
380            momentum_y = quantities[3]
381            depth = w - z
382             
383            if w == NAN or z == NAN or momentum_x == NAN:
384                velocity_x = NAN
385            else:
386                if depth > 1.e-30: # use epsilon
387                    velocity_x = momentum_x / depth  #Absolute velocity
388                else:
389                    velocity_x = 0
390            if w == NAN or z == NAN or momentum_y == NAN:
391                velocity_y = NAN
392            else:
393                if depth > 1.e-30: # use epsilon
394                    velocity_y = momentum_y / depth  #Absolute velocity
395                else:
396                    velocity_y = 0
397            depths.append(depth)
398            velocity_xs.append(velocity_x)
399            velocity_ys.append(velocity_y)
400            if stage_file is not None:
401                stages.append(w)
402        depth_writer.writerow(depths)
403        velocity_x_writer.writerow(velocity_xs)
404        velocity_y_writer.writerow(velocity_ys)
405        if stage_file is not None:
406            stage_writer.writerow(stages)       
407
408
409class Interpolation_function:
410    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
411    which is interpolated from time series defined at vertices of
412    triangular mesh (such as those stored in sww files)
413
414    Let m be the number of vertices, n the number of triangles
415    and p the number of timesteps.
416
417    Mandatory input
418        time:               px1 array of monotonously increasing times (Float)
419        quantities:         Dictionary of arrays or 1 array (Float)
420                            The arrays must either have dimensions pxm or mx1.
421                            The resulting function will be time dependent in
422                            the former case while it will be constan with
423                            respect to time in the latter case.
424       
425    Optional input:
426        quantity_names:     List of keys into the quantities dictionary
427        vertex_coordinates: mx2 array of coordinates (Float)
428        triangles:          nx3 array of indices into vertex_coordinates (Int)
429        interpolation_points: Nx2 array of coordinates to be interpolated to
430        verbose:            Level of reporting
431   
432   
433    The quantities returned by the callable object are specified by
434    the list quantities which must contain the names of the
435    quantities to be returned and also reflect the order, e.g. for
436    the shallow water wave equation, on would have
437    quantities = ['stage', 'xmomentum', 'ymomentum']
438
439    The parameter interpolation_points decides at which points interpolated
440    quantities are to be computed whenever object is called.
441    If None, return average value
442
443    FIXME (Ole): Need to allow vertex coordinates and interpolation points to be
444    geospatial data objects
445
446    Time assumed to be relative to starttime   
447    """
448 
449   
450    def __init__(self,
451                 time,
452                 quantities,
453                 quantity_names=None, 
454                 vertex_coordinates=None,
455                 triangles=None,
456                 interpolation_points=None,
457                 time_thinning=1,
458                 verbose=False):
459        """Initialise object and build spatial interpolation if required
460
461        Time_thinning_number controls how many timesteps to use. Only timesteps with
462        index%time_thinning_number == 0 will used, or in other words a value of 3, say,
463        will cause the algorithm to use every third time step.
464        """
465
466        from Numeric import array, zeros, Float, alltrue, concatenate,\
467             reshape, ArrayType
468
469
470        from anuga.config import time_format
471        import types
472
473
474        # Check temporal info
475        time = ensure_numeric(time)       
476        msg = 'Time must be a monotonuosly '
477        msg += 'increasing sequence %s' %time
478        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
479
480
481        # Check if quantities is a single array only
482        if type(quantities) != types.DictType:
483            quantities = ensure_numeric(quantities)
484            quantity_names = ['Attribute']
485
486            #Make it a dictionary
487            quantities = {quantity_names[0]: quantities}
488
489
490        # Use keys if no names are specified
491        if quantity_names is None:
492            quantity_names = quantities.keys()
493
494
495        # Check spatial info
496        if vertex_coordinates is None:
497            self.spatial = False
498        else:   
499            vertex_coordinates = ensure_absolute(vertex_coordinates)
500
501            assert triangles is not None, 'Triangles array must be specified'
502            triangles = ensure_numeric(triangles)
503            self.spatial = True           
504
505        # Thin timesteps if needed
506        # Note array() is used to make the thinned arrays contiguous in memory
507        self.time = array(time[::time_thinning])         
508        for name in quantity_names:
509            if len(quantities[name].shape) == 2:
510                quantities[name] = array(quantities[name][::time_thinning,:])
511             
512        # Save for use with statistics
513        self.quantities_range = {}
514        for name in quantity_names:
515            q = quantities[name][:].flat
516            self.quantities_range[name] = [min(q), max(q)]
517       
518        self.quantity_names = quantity_names       
519        self.vertex_coordinates = vertex_coordinates
520        self.interpolation_points = interpolation_points
521       
522
523        self.index = 0    # Initial time index
524        self.precomputed_values = {}
525       
526           
527        # Precomputed spatial interpolation if requested
528        if interpolation_points is not None:
529            if self.spatial is False:
530                raise 'Triangles and vertex_coordinates must be specified'
531           
532            try:
533                self.interpolation_points = interpolation_points = ensure_numeric(interpolation_points)
534            except:
535                msg = 'Interpolation points must be an N x 2 Numeric array '+\
536                      'or a list of points\n'
537                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
538                                      '...')
539                raise msg
540
541
542            # Check that all interpolation points fall within
543            # mesh boundary as defined by triangles and vertex_coordinates.
544            from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
545            from anuga.utilities.polygon import outside_polygon           
546
547            # Create temporary mesh object from mesh info passed
548            # into this function.
549            mesh = Mesh(vertex_coordinates, triangles)
550            mesh_boundary_polygon = mesh.get_boundary_polygon()
551
552           
553            indices = outside_polygon(interpolation_points,
554                                      mesh_boundary_polygon)
555
556            # Record result
557            #self.mesh_boundary_polygon = mesh_boundary_polygon
558            self.indices_outside_mesh = indices
559
560            # Report
561            if len(indices) > 0:
562                msg = 'Interpolation points in Interpolation function fall ' 
563                msg += 'outside specified mesh. '
564                msg += 'Offending points:\n'
565                out_interp_pts = []
566                for i in indices:
567                    msg += '%d: %s\n' %(i, interpolation_points[i])
568                    out_interp_pts.append(ensure_numeric(interpolation_points[i]))
569
570                if verbose is True:
571                    import sys
572                    if sys.platform == 'win32':
573                        from anuga.utilities.polygon import plot_polygons
574                        #out_interp_pts = take(interpolation_points,[indices])
575                        title = 'Interpolation points fall outside specified mesh'
576                        plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts],
577                                      ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose)
578
579                # Joaquim Luis suggested this as an Exception, so
580                # that the user can now what the problem is rather than
581                # looking for NaN's. However, NANs are handy as they can
582                # be ignored leaving good points for continued processing.
583                if verbose:
584                    print msg
585                #raise Exception(msg)
586
587            # Plot boundary and interpolation points
588            if verbose is True:
589                import sys
590                if sys.platform == 'win32':
591                    from anuga.utilities.polygon import plot_polygons
592                    title = 'Interpolation function: Polygon and interpolation points'
593                    plot_polygons([mesh_boundary_polygon,interpolation_points],
594                                         ['line','point'],figname='points_boundary',label=title,verbose=verbose)
595
596            m = len(self.interpolation_points)
597            p = len(self.time)
598           
599            for name in quantity_names:
600                self.precomputed_values[name] = zeros((p, m), Float)
601
602            # Build interpolator
603            if verbose:
604                msg = 'Building interpolation matrix from source mesh '
605                msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
606                                                       triangles.shape[0])
607                print msg
608
609               
610            interpol = Interpolate(vertex_coordinates,
611                                   triangles,
612                                   verbose=verbose)
613
614            if verbose:
615                print 'Interpolating (%d interpolation points, %d timesteps).'\
616                      %(self.interpolation_points.shape[0], self.time.shape[0]),
617           
618                if time_thinning > 1:
619                    print 'Timesteps were thinned by a factor of %d' %time_thinning
620                else:
621                    print
622
623            for i, t in enumerate(self.time):
624                # Interpolate quantities at this timestep
625                if verbose and i%((p+10)/10)==0:
626                    print ' time step %d of %d' %(i, p)
627                   
628                for name in quantity_names:
629                    if len(quantities[name].shape) == 2:
630                        Q = quantities[name][i,:] # Quantities at timestep i
631                    else:
632                        Q = quantities[name][:]   # No time dependency
633                       
634                    # Interpolate   
635                    result = interpol.interpolate(Q,
636                                                  point_coordinates=\
637                                                  self.interpolation_points)
638
639                    #assert len(result), len(interpolation_points)
640                    self.precomputed_values[name][i, :] = result
641
642                   
643            # Report
644            if verbose:
645                print self.statistics()
646                #self.print_statistics()
647           
648        else:
649            # Store quantitites as is
650            for name in quantity_names:
651                self.precomputed_values[name] = quantities[name]
652
653    def __repr__(self):
654        # return 'Interpolation function (spatio-temporal)'
655        return self.statistics()
656 
657    def __call__(self, t, point_id=None, x=None, y=None):
658        """Evaluate f(t) or f(t, point_id)
659       
660        Inputs:
661          t: time - Model time. Must lie within existing timesteps
662          point_id: index of one of the preprocessed points.
663     
664                   
665          If spatial info is present and all of point_id
666          are None an exception is raised
667                   
668          If no spatial info is present, point_id arguments are ignored
669          making f a function of time only.
670
671
672          FIXME: f(t, x, y) x, y could overrided location, point_id ignored 
673          FIXME: point_id could also be a slice
674          FIXME: What if x and y are vectors?
675          FIXME: What about f(x,y) without t?
676        """
677
678        from math import pi, cos, sin, sqrt
679        from Numeric import zeros, Float
680        from anuga.abstract_2d_finite_volumes.util import mean       
681
682        if self.spatial is True:
683            if point_id is None:
684                if x is None or y is None:
685                    msg = 'Either point_id or x and y must be specified'
686                    raise Exception(msg)
687            else:
688                if self.interpolation_points is None:
689                    msg = 'Interpolation_function must be instantiated ' +\
690                          'with a list of interpolation points before parameter ' +\
691                          'point_id can be used'
692                    raise Exception(msg)
693
694        msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
695        msg += ' does not match model time: %.16f\n' %t
696        if t < self.time[0]: raise Exception(msg)
697        if t > self.time[-1]: raise Exception(msg)
698
699        oldindex = self.index #Time index
700
701        # Find current time slot
702        while t > self.time[self.index]: self.index += 1
703        while t < self.time[self.index]: self.index -= 1
704
705        if t == self.time[self.index]:
706            # Protect against case where t == T[-1] (last time)
707            #  - also works in general when t == T[i]
708            ratio = 0
709        else:
710            # t is now between index and index+1
711            ratio = (t - self.time[self.index])/\
712                    (self.time[self.index+1] - self.time[self.index])
713
714        # Compute interpolated values
715        q = zeros(len(self.quantity_names), Float)
716        # print "self.precomputed_values", self.precomputed_values
717        for i, name in enumerate(self.quantity_names):
718            Q = self.precomputed_values[name]
719
720            if self.spatial is False:
721                # If there is no spatial info               
722                assert len(Q.shape) == 1
723
724                Q0 = Q[self.index]
725                if ratio > 0: Q1 = Q[self.index+1]
726
727            else:
728                if x is not None and y is not None:
729                    # Interpolate to x, y
730                   
731                    raise 'x,y interpolation not yet implemented'
732                else:
733                    # Use precomputed point
734                    Q0 = Q[self.index, point_id]
735                    if ratio > 0:
736                        Q1 = Q[self.index+1, point_id]
737
738            # Linear temporal interpolation   
739            if ratio > 0:
740                if Q0 == NAN and Q1 == NAN:
741                    q[i]  = Q0
742                else:
743                    q[i] = Q0 + ratio*(Q1 - Q0)
744            else:
745                q[i] = Q0
746
747
748        # Return vector of interpolated values
749        # if len(q) == 1:
750        #     return q[0]
751        # else:
752        #     return q
753
754
755        # Return vector of interpolated values
756        # FIXME:
757        if self.spatial is True:
758            return q
759        else:
760            # Replicate q according to x and y
761            # This is e.g used for Wind_stress
762            if x is None or y is None: 
763                return q
764            else:
765                try:
766                    N = len(x)
767                except:
768                    return q
769                else:
770                    from Numeric import ones, Float
771                    # x is a vector - Create one constant column for each value
772                    N = len(x)
773                    assert len(y) == N, 'x and y must have same length'
774                    res = []
775                    for col in q:
776                        res.append(col*ones(N, Float))
777                       
778                return res
779
780
781    def get_time(self):
782        """Return model time as a vector of timesteps
783        """
784        return self.time
785
786
787    def statistics(self):
788        """Output statistics about interpolation_function
789        """
790       
791        vertex_coordinates = self.vertex_coordinates
792        interpolation_points = self.interpolation_points               
793        quantity_names = self.quantity_names
794        #quantities = self.quantities
795        precomputed_values = self.precomputed_values                 
796               
797        x = vertex_coordinates[:,0]
798        y = vertex_coordinates[:,1]               
799
800        str =  '------------------------------------------------\n'
801        str += 'Interpolation_function (spatio-temporal) statistics:\n'
802        str += '  Extent:\n'
803        str += '    x in [%f, %f], len(x) == %d\n'\
804               %(min(x), max(x), len(x))
805        str += '    y in [%f, %f], len(y) == %d\n'\
806               %(min(y), max(y), len(y))
807        str += '    t in [%f, %f], len(t) == %d\n'\
808               %(min(self.time), max(self.time), len(self.time))
809        str += '  Quantities:\n'
810        for name in quantity_names:
811            minq, maxq = self.quantities_range[name]
812            str += '    %s in [%f, %f]\n' %(name, minq, maxq)           
813            #q = quantities[name][:].flat
814            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
815
816        if interpolation_points is not None:   
817            str += '  Interpolation points (xi, eta):'\
818                   ' number of points == %d\n' %interpolation_points.shape[0]
819            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
820                                            max(interpolation_points[:,0]))
821            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
822                                             max(interpolation_points[:,1]))
823            str += '  Interpolated quantities (over all timesteps):\n'
824       
825            for name in quantity_names:
826                q = precomputed_values[name][:].flat
827                str += '    %s at interpolation points in [%f, %f]\n'\
828                       %(name, min(q), max(q))
829        str += '------------------------------------------------\n'
830
831        return str
832
833
834def interpolate_sww(sww_file, time, interpolation_points,
835                    quantity_names=None, verbose=False):
836    """
837    obsolete.
838    use file_function in utils
839    """
840    #open sww file
841    x, y, volumes, time, quantities = read_sww(sww_file)
842    print "x",x
843    print "y",y
844   
845    print "time", time
846    print "quantities", quantities
847
848    #Add the x and y together
849    vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
850
851    #Will return the quantity values at the specified times and locations
852    interp = Interpolation_interface(time,
853                                     quantities,
854                                     quantity_names=quantity_names, 
855                                     vertex_coordinates=vertex_coordinates,
856                                     triangles=volumes,
857                                     interpolation_points=interpolation_points,
858                                     verbose=verbose)
859
860
861def read_sww(file_name):
862    """
863    obsolete - Nothing should be calling this
864   
865    Read in an sww file.
866   
867    Input;
868    file_name - the sww file
869   
870    Output;
871    x - Vector of x values
872    y - Vector of y values
873    z - Vector of bed elevation
874    volumes - Array.  Each row has 3 values, representing
875    the vertices that define the volume
876    time - Vector of the times where there is stage information
877    stage - array with respect to time and vertices (x,y)
878    """
879
880    msg = 'Function read_sww in interpolat.py is obsolete'
881    raise Exception, msg
882   
883    #FIXME Have this reader as part of data_manager?
884   
885    from Scientific.IO.NetCDF import NetCDFFile     
886    import tempfile
887    import sys
888    import os
889   
890    #Check contents
891    #Get NetCDF
892   
893    # see if the file is there.  Throw a QUIET IO error if it isn't
894    # I don't think I could implement the above
895   
896    #throws prints to screen if file not present
897    junk = tempfile.mktemp(".txt")
898    fd = open(junk,'w')
899    stdout = sys.stdout
900    sys.stdout = fd
901    fid = NetCDFFile(file_name, 'r') 
902    sys.stdout = stdout
903    fd.close()
904    #clean up
905    os.remove(junk)     
906     
907    # Get the variables
908    x = fid.variables['x'][:]
909    y = fid.variables['y'][:]
910    volumes = fid.variables['volumes'][:] 
911    time = fid.variables['time'][:]
912
913    keys = fid.variables.keys()
914    keys.remove("x")
915    keys.remove("y")
916    keys.remove("volumes")
917    keys.remove("time")
918     #Turn NetCDF objects into Numeric arrays
919    quantities = {}
920    for name in keys:
921        quantities[name] = fid.variables[name][:]
922           
923    fid.close()
924    return x, y, volumes, time, quantities
925
926
927#-------------------------------------------------------------
928if __name__ == "__main__":
929    names = ["x","y"]
930    someiterable = [[1,2],[3,4]]
931    csvwriter = writer(file("some.csv", "wb"))
932    csvwriter.writerow(names)
933    for row in someiterable:
934        csvwriter.writerow(row)
Note: See TracBrowser for help on using the repository browser.