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

Last change on this file since 3900 was 3900, checked in by ole, 17 years ago

Added time_thinning to Interpolation_function and improved diagnostics.

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