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

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

reduce memory use in quantity.set_value. fit_to_mesh can now use an existing mesh instance, which it does in quantity.set_value.

File size: 32.3 KB
Line 
1"""Least squares interpolation.
2
3   Implements a least-squares interpolation.
4   Putting mesh data onto points.
5
6   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
7   Geoscience Australia, 2004.
8
9DESIGN ISSUES
10* what variables should be global?
11- if there are no global vars functions can be moved around alot easier
12
13* The public interface to Interpolate
14__init__
15interpolate
16interpolate_block
17
18"""
19
20import time
21import os
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    """
299    No test for this yet.
300    Note, this has no time the input data has no time dimension.  Which is
301    different from most of the data we interpolate, eg sww info.
302     
303        Output:
304          Interpolated values at inputted points.
305    """
306    interp = Interpolate(vertices,
307                         triangles, 
308                         max_vertices_per_cell=max_points_per_cell)
309           
310    calc = interp.interpolate(vertex_attributes
311                              ,points
312                              ,start_blocking_len=start_blocking_len)
313def interpolate_sww2csv(sww_file,
314                        points,
315                        depth_file,
316                        velocity_x_file,
317                        velocity_y_file,
318                        #quantities = ['depth', 'velocity'],
319                        verbose=True,
320                        use_cache = True):
321    """
322    Interpolate the quantities at a given set of locations, given
323    an sww file.
324    The results are written to a csv file.
325
326    In the future let points be a points file.
327    And the user choose the quantities.
328
329    This is currently quite specific.
330    If it need to be more general, chagne things.
331
332    This is really returning speed, not velocity.
333    """
334    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
335    #print "points",points
336    points = ensure_absolute(points)
337    point_count = len(points)
338    callable_sww = file_function(sww_file,
339                                 quantities=quantities,
340                                 interpolation_points=points,
341                                 verbose=verbose,
342                                 use_cache=use_cache)
343   
344    depth_writer = writer(file(depth_file, "wb"))
345    velocity_x_writer = writer(file(velocity_x_file, "wb"))
346    velocity_y_writer = writer(file(velocity_y_file, "wb"))
347    # Write heading
348    heading = [str(x[0])+ ':' + str(x[1]) for x in points]
349    heading.insert(0, "time")
350    depth_writer.writerow(heading)
351    velocity_x_writer.writerow(heading)
352    velocity_y_writer.writerow(heading)
353   
354    for time in callable_sww.get_time():
355        depths = [time]
356        velocity_xs = [time]
357        velocity_ys = [time]
358        for point_i, point in enumerate(points):
359            quantities = callable_sww(time,point_i)
360            #print "quantities", quantities
361           
362            w = quantities[0]
363            z = quantities[1]
364            momentum_x = quantities[2]
365            momentum_y = quantities[3]
366            depth = w - z
367             
368            if w == NAN or z == NAN or momentum_x == NAN:
369                velocity_x = NAN
370            else:
371                if depth > 1.e-30: # use epsilon
372                    velocity_x = momentum_x / depth  #Absolute velocity
373                else:
374                    velocity_x = 0
375            if w == NAN or z == NAN or momentum_y == NAN:
376                velocity_y = NAN
377            else:
378                if depth > 1.e-30: # use epsilon
379                    velocity_y = momentum_y / depth  #Absolute velocity
380                else:
381                    velocity_y = 0
382            depths.append(depth)
383            velocity_xs.append(velocity_x)
384            velocity_ys.append(velocity_y)
385        depth_writer.writerow(depths)
386        velocity_x_writer.writerow(velocity_xs)
387        velocity_y_writer.writerow(velocity_ys)
388
389
390class Interpolation_function:
391    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
392    which is interpolated from time series defined at vertices of
393    triangular mesh (such as those stored in sww files)
394
395    Let m be the number of vertices, n the number of triangles
396    and p the number of timesteps.
397
398    Mandatory input
399        time:               px1 array of monotonously increasing times (Float)
400        quantities:         Dictionary of arrays or 1 array (Float)
401                            The arrays must either have dimensions pxm or mx1.
402                            The resulting function will be time dependent in
403                            the former case while it will be constan with
404                            respect to time in the latter case.
405       
406    Optional input:
407        quantity_names:     List of keys into the quantities dictionary
408        vertex_coordinates: mx2 array of coordinates (Float)
409        triangles:          nx3 array of indices into vertex_coordinates (Int)
410        interpolation_points: Nx2 array of coordinates to be interpolated to
411        verbose:            Level of reporting
412   
413   
414    The quantities returned by the callable object are specified by
415    the list quantities which must contain the names of the
416    quantities to be returned and also reflect the order, e.g. for
417    the shallow water wave equation, on would have
418    quantities = ['stage', 'xmomentum', 'ymomentum']
419
420    The parameter interpolation_points decides at which points interpolated
421    quantities are to be computed whenever object is called.
422    If None, return average value
423
424    FIXME (Ole): Need to allow vertex coordinates and interpolation points to be
425    geospatial data objects
426
427    Time assumed to be relative to starttime   
428    """
429 
430   
431    def __init__(self,
432                 time,
433                 quantities,
434                 quantity_names=None, 
435                 vertex_coordinates=None,
436                 triangles=None,
437                 interpolation_points=None,
438                 time_thinning=1,
439                 verbose=False):
440        """Initialise object and build spatial interpolation if required
441
442        Time_thinning_number controls how many timesteps to use. Only timesteps with
443        index%time_thinning_number == 0 will used, or in other words a value of 3, say,
444        will cause the algorithm to use every third time step.
445        """
446
447        from Numeric import array, zeros, Float, alltrue, concatenate,\
448             reshape, ArrayType
449
450
451        from anuga.config import time_format
452        import types
453
454
455        # Check temporal info
456        time = ensure_numeric(time)       
457        msg = 'Time must be a monotonuosly '
458        msg += 'increasing sequence %s' %time
459        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
460
461
462        # Check if quantities is a single array only
463        if type(quantities) != types.DictType:
464            quantities = ensure_numeric(quantities)
465            quantity_names = ['Attribute']
466
467            #Make it a dictionary
468            quantities = {quantity_names[0]: quantities}
469
470
471        # Use keys if no names are specified
472        if quantity_names is None:
473            quantity_names = quantities.keys()
474
475
476        # Check spatial info
477        if vertex_coordinates is None:
478            self.spatial = False
479        else:   
480            vertex_coordinates = ensure_absolute(vertex_coordinates)
481
482            assert triangles is not None, 'Triangles array must be specified'
483            triangles = ensure_numeric(triangles)
484            self.spatial = True           
485
486        # Thin timesteps if needed
487        # Note array() is used to make the thinned arrays contiguous in memory
488        self.time = array(time[::time_thinning])         
489        for name in quantity_names:
490            if len(quantities[name].shape) == 2:
491                quantities[name] = array(quantities[name][::time_thinning,:])
492             
493        # Save for use with statistics
494        self.quantities_range = {}
495        for name in quantity_names:
496            q = quantities[name][:].flat
497            self.quantities_range[name] = [min(q), max(q)]
498       
499        self.quantity_names = quantity_names       
500        self.vertex_coordinates = vertex_coordinates
501        self.interpolation_points = interpolation_points
502       
503
504        self.index = 0    # Initial time index
505        self.precomputed_values = {}
506       
507           
508        # Precomputed spatial interpolation if requested
509        if interpolation_points is not None:
510            if self.spatial is False:
511                raise 'Triangles and vertex_coordinates must be specified'
512           
513            try:
514                self.interpolation_points = ensure_numeric(interpolation_points)
515            except:
516                msg = 'Interpolation points must be an N x 2 Numeric array '+\
517                      'or a list of points\n'
518                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
519                                      '...')
520                raise msg
521
522
523            m = len(self.interpolation_points)
524            p = len(self.time)
525           
526            for name in quantity_names:
527                self.precomputed_values[name] = zeros((p, m), Float)
528
529            # Build interpolator
530            if verbose:
531                msg = 'Building interpolation matrix from source mesh '
532                msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
533                                                       triangles.shape[0])
534                print msg
535
536               
537            interpol = Interpolate(vertex_coordinates,
538                                   triangles,
539                                   verbose=verbose)
540
541            if verbose:
542                print 'Interpolating (%d interpolation points, %d timesteps).'\
543                      %(self.interpolation_points.shape[0], self.time.shape[0]),
544           
545                if time_thinning > 1:
546                    print 'Timesteps were thinned by a factor of %d' %time_thinning
547                else:
548                    print
549
550            for i, t in enumerate(self.time):
551                # Interpolate quantities at this timestep
552                if verbose and i%((p+10)/10)==0:
553                    print ' time step %d of %d' %(i, p)
554                   
555                for name in quantity_names:
556                    if len(quantities[name].shape) == 2:
557                        Q = quantities[name][i,:] # Quantities at timestep i
558                    else:
559                        Q = quantities[name][:]   # No time dependency
560                       
561                    # Interpolate   
562                    result = interpol.interpolate(Q,
563                                                  point_coordinates=\
564                                                  self.interpolation_points)
565
566                    #assert len(result), len(interpolation_points)
567                    self.precomputed_values[name][i, :] = result
568
569                   
570            # Report
571            if verbose:
572                print self.statistics()
573                #self.print_statistics()
574           
575        else:
576            # Store quantitites as is
577            for name in quantity_names:
578                self.precomputed_values[name] = quantities[name]
579
580    def __repr__(self):
581        # return 'Interpolation function (spatio-temporal)'
582        return self.statistics()
583 
584    def __call__(self, t, point_id=None, x=None, y=None):
585        """Evaluate f(t) or f(t, point_id)
586       
587        Inputs:
588          t: time - Model time. Must lie within existing timesteps
589          point_id: index of one of the preprocessed points.
590     
591                   
592          If spatial info is present and all of point_id
593          are None an exception is raised
594                   
595          If no spatial info is present, point_id arguments are ignored
596          making f a function of time only.
597
598
599          FIXME: f(t, x, y) x, y could overrided location, point_id ignored 
600          FIXME: point_id could also be a slice
601          FIXME: What if x and y are vectors?
602          FIXME: What about f(x,y) without t?
603        """
604
605        from math import pi, cos, sin, sqrt
606        from Numeric import zeros, Float
607        from anuga.abstract_2d_finite_volumes.util import mean       
608
609        if self.spatial is True:
610            if point_id is None:
611                if x is None or y is None:
612                    msg = 'Either point_id or x and y must be specified'
613                    raise Exception(msg)
614            else:
615                if self.interpolation_points is None:
616                    msg = 'Interpolation_function must be instantiated ' +\
617                          'with a list of interpolation points before parameter ' +\
618                          'point_id can be used'
619                    raise Exception(msg)
620
621        msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
622        msg += ' does not match model time: %.16f\n' %t
623        if t < self.time[0]: raise Exception(msg)
624        if t > self.time[-1]: raise Exception(msg)
625
626        oldindex = self.index #Time index
627
628        #Find current time slot
629        while t > self.time[self.index]: self.index += 1
630        while t < self.time[self.index]: self.index -= 1
631
632        if t == self.time[self.index]:
633            #Protect against case where t == T[-1] (last time)
634            # - also works in general when t == T[i]
635            ratio = 0
636        else:
637            #t is now between index and index+1
638            ratio = (t - self.time[self.index])/\
639                    (self.time[self.index+1] - self.time[self.index])
640
641        #Compute interpolated values
642        q = zeros(len(self.quantity_names), Float)
643        #print "self.precomputed_values", self.precomputed_values
644        for i, name in enumerate(self.quantity_names):
645            Q = self.precomputed_values[name]
646
647            if self.spatial is False:
648                #If there is no spatial info               
649                assert len(Q.shape) == 1
650
651                Q0 = Q[self.index]
652                if ratio > 0: Q1 = Q[self.index+1]
653
654            else:
655                if x is not None and y is not None:
656                    #Interpolate to x, y
657                   
658                    raise 'x,y interpolation not yet implemented'
659                else:
660                    #Use precomputed point
661                    Q0 = Q[self.index, point_id]
662                    if ratio > 0:
663                        Q1 = Q[self.index+1, point_id]
664
665            #Linear temporal interpolation   
666            if ratio > 0:
667                if Q0 == NAN and Q1 == NAN:
668                    q[i]  = Q0
669                else:
670                    q[i] = Q0 + ratio*(Q1 - Q0)
671            else:
672                q[i] = Q0
673
674
675        #Return vector of interpolated values
676        #if len(q) == 1:
677        #    return q[0]
678        #else:
679        #    return q
680
681
682        #Return vector of interpolated values
683        #FIXME:
684        if self.spatial is True:
685            return q
686        else:
687            #Replicate q according to x and y
688            #This is e.g used for Wind_stress
689            if x is None or y is None: 
690                return q
691            else:
692                try:
693                    N = len(x)
694                except:
695                    return q
696                else:
697                    from Numeric import ones, Float
698                    #x is a vector - Create one constant column for each value
699                    N = len(x)
700                    assert len(y) == N, 'x and y must have same length'
701                    res = []
702                    for col in q:
703                        res.append(col*ones(N, Float))
704                       
705                return res
706
707
708    def get_time(self):
709        """Return model time as a vector of timesteps
710        """
711        return self.time
712
713
714    def statistics(self):
715        """Output statistics about interpolation_function
716        """
717       
718        vertex_coordinates = self.vertex_coordinates
719        interpolation_points = self.interpolation_points               
720        quantity_names = self.quantity_names
721        #quantities = self.quantities
722        precomputed_values = self.precomputed_values                 
723               
724        x = vertex_coordinates[:,0]
725        y = vertex_coordinates[:,1]               
726
727        str =  '------------------------------------------------\n'
728        str += 'Interpolation_function (spatio-temporal) statistics:\n'
729        str += '  Extent:\n'
730        str += '    x in [%f, %f], len(x) == %d\n'\
731               %(min(x), max(x), len(x))
732        str += '    y in [%f, %f], len(y) == %d\n'\
733               %(min(y), max(y), len(y))
734        str += '    t in [%f, %f], len(t) == %d\n'\
735               %(min(self.time), max(self.time), len(self.time))
736        str += '  Quantities:\n'
737        for name in quantity_names:
738            minq, maxq = self.quantities_range[name]
739            str += '    %s in [%f, %f]\n' %(name, minq, maxq)           
740            #q = quantities[name][:].flat
741            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
742
743        if interpolation_points is not None:   
744            str += '  Interpolation points (xi, eta):'\
745                   ' number of points == %d\n' %interpolation_points.shape[0]
746            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
747                                            max(interpolation_points[:,0]))
748            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
749                                             max(interpolation_points[:,1]))
750            str += '  Interpolated quantities (over all timesteps):\n'
751       
752            for name in quantity_names:
753                q = precomputed_values[name][:].flat
754                str += '    %s at interpolation points in [%f, %f]\n'\
755                       %(name, min(q), max(q))
756        str += '------------------------------------------------\n'
757
758        return str
759
760
761def interpolate_sww(sww_file, time, interpolation_points,
762                    quantity_names=None, verbose=False):
763    """
764    obsolete.
765    use file_function in utils
766    """
767    #open sww file
768    x, y, volumes, time, quantities = read_sww(sww_file)
769    print "x",x
770    print "y",y
771   
772    print "time", time
773    print "quantities", quantities
774
775    #Add the x and y together
776    vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
777
778    #Will return the quantity values at the specified times and locations
779    interp = Interpolation_interface(time,
780                                     quantities,
781                                     quantity_names=quantity_names, 
782                                     vertex_coordinates=vertex_coordinates,
783                                     triangles=volumes,
784                                     interpolation_points=interpolation_points,
785                                     verbose=verbose)
786
787
788def read_sww(file_name):
789    """
790    obsolete - Nothing should be calling this
791   
792    Read in an sww file.
793   
794    Input;
795    file_name - the sww file
796   
797    Output;
798    x - Vector of x values
799    y - Vector of y values
800    z - Vector of bed elevation
801    volumes - Array.  Each row has 3 values, representing
802    the vertices that define the volume
803    time - Vector of the times where there is stage information
804    stage - array with respect to time and vertices (x,y)
805    """
806
807    msg = 'Function read_sww in interpolat.py is obsolete'
808    raise Exception, msg
809   
810    #FIXME Have this reader as part of data_manager?
811   
812    from Scientific.IO.NetCDF import NetCDFFile     
813    import tempfile
814    import sys
815    import os
816   
817    #Check contents
818    #Get NetCDF
819   
820    # see if the file is there.  Throw a QUIET IO error if it isn't
821    # I don't think I could implement the above
822   
823    #throws prints to screen if file not present
824    junk = tempfile.mktemp(".txt")
825    fd = open(junk,'w')
826    stdout = sys.stdout
827    sys.stdout = fd
828    fid = NetCDFFile(file_name, 'r') 
829    sys.stdout = stdout
830    fd.close()
831    #clean up
832    os.remove(junk)     
833     
834    # Get the variables
835    x = fid.variables['x'][:]
836    y = fid.variables['y'][:]
837    volumes = fid.variables['volumes'][:] 
838    time = fid.variables['time'][:]
839
840    keys = fid.variables.keys()
841    keys.remove("x")
842    keys.remove("y")
843    keys.remove("volumes")
844    keys.remove("time")
845     #Turn NetCDF objects into Numeric arrays
846    quantities = {}
847    for name in keys:
848        quantities[name] = fid.variables[name][:]
849           
850    fid.close()
851    return x, y, volumes, time, quantities
852
853
854#-------------------------------------------------------------
855if __name__ == "__main__":
856    names = ["x","y"]
857    someiterable = [[1,2],[3,4]]
858    csvwriter = writer(file("some.csv", "wb"))
859    csvwriter.writerow(names)
860    for row in someiterable:
861        csvwriter.writerow(row)
Note: See TracBrowser for help on using the repository browser.