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

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

Added functionality for getting arbitrary interpolated values in Quantity as well as calculating inundation height and location. This work was done at SUT during the last week of September 2006.

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