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

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

Comments

File size: 30.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 
409   
410    def __init__(self,
411                 time,
412                 quantities,
413                 quantity_names=None, 
414                 vertex_coordinates=None,
415                 triangles=None,
416                 interpolation_points=None,
417                 verbose=False):
418        """Initialise object and build spatial interpolation if required
419        """
420
421        from Numeric import array, zeros, Float, alltrue, concatenate,\
422             reshape, ArrayType
423
424
425        from anuga.config import time_format
426        import types
427
428
429        # Check temporal info
430        time = ensure_numeric(time)       
431        msg = 'Time must be a monotonuosly '
432        msg += 'increasing sequence %s' %time
433        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
434
435
436        # Check if quantities is a single array only
437        if type(quantities) != types.DictType:
438            quantities = ensure_numeric(quantities)
439            quantity_names = ['Attribute']
440
441            #Make it a dictionary
442            quantities = {quantity_names[0]: quantities}
443
444
445        # Use keys if no names are specified
446        if quantity_names is None:
447            quantity_names = quantities.keys()
448
449
450        # Check spatial info
451        if vertex_coordinates is None:
452            self.spatial = False
453        else:   
454            vertex_coordinates = ensure_numeric(vertex_coordinates)
455
456            assert triangles is not None, 'Triangles array must be specified'
457            triangles = ensure_numeric(triangles)
458            self.spatial = True           
459
460             
461        # Save for use with statistics
462        self.quantities_range = {}
463        for name in quantity_names:
464            q = quantities[name][:].flat
465            self.quantities_range[name] = [min(q), max(q)]
466       
467        self.quantity_names = quantity_names       
468        self.vertex_coordinates = vertex_coordinates
469        self.interpolation_points = interpolation_points
470        self.time = time[:]  # Time assumed to be relative to starttime
471        self.index = 0    # Initial time index
472        self.precomputed_values = {}
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            if verbose: print 'Building interpolation matrix'           
498            interpol = Interpolate(vertex_coordinates,
499                                   triangles,
500                                   #point_coordinates = \
501                                   #self.interpolation_points,
502                                   #alpha = 0,
503                                   verbose=verbose)
504
505            if verbose: print 'Interpolate'
506            for i, t in enumerate(self.time):
507                # Interpolate quantities at this timestep
508                if verbose and i%((p+10)/10)==0:
509                    print ' time step %d of %d' %(i, p)
510                   
511                for name in quantity_names:
512                    if len(quantities[name].shape) == 2:
513                        Q = quantities[name][i,:] # Quantities at timestep i
514                    else:
515                        Q = quantities[name][:]   # No time dependency
516                       
517                    # Interpolate   
518                    result = interpol.interpolate(Q,
519                                                  point_coordinates=\
520                                                  self.interpolation_points)
521                    self.precomputed_values[name][i, :] = result
522
523                   
524            # Report
525            if verbose:
526                print self.statistics()
527                #self.print_statistics()
528           
529        else:
530            # Store quantitites as is
531            for name in quantity_names:
532                self.precomputed_values[name] = quantities[name]
533
534    def __repr__(self):
535        # return 'Interpolation function (spatio-temporal)'
536        return self.statistics()
537 
538    def __call__(self, t, point_id=None, x=None, y=None):
539        """Evaluate f(t), f(t, point_id) or f(t, x, y)
540
541        Inputs:
542          t: time - Model time. Must lie within existing timesteps
543          point_id: index of one of the preprocessed points.
544          x, y:     Overrides location, point_id ignored
545         
546          If spatial info is present and all of x,y,point_id
547          are None an exception is raised
548                   
549          If no spatial info is present, point_id and x,y arguments are ignored
550          making f a function of time only.
551
552         
553          FIXME: point_id could also be a slice
554          FIXME: What if x and y are vectors?
555          FIXME: What about f(x,y) without t?
556        """
557
558        from math import pi, cos, sin, sqrt
559        from Numeric import zeros, Float
560        from anuga.abstract_2d_finite_volumes.util import mean       
561
562        if self.spatial is True:
563            if point_id is None:
564                if x is None or y is None:
565                    msg = 'Either point_id or x and y must be specified'
566                    raise Exception(msg)
567            else:
568                if self.interpolation_points is None:
569                    msg = 'Interpolation_function must be instantiated ' +\
570                          'with a list of interpolation points before parameter ' +\
571                          'point_id can be used'
572                    raise Exception(msg)
573
574        msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
575        msg += ' does not match model time: %.16f\n' %t
576        if t < self.time[0]: raise Exception(msg)
577        if t > self.time[-1]: raise Exception(msg)
578
579        oldindex = self.index #Time index
580
581        #Find current time slot
582        while t > self.time[self.index]: self.index += 1
583        while t < self.time[self.index]: self.index -= 1
584
585        if t == self.time[self.index]:
586            #Protect against case where t == T[-1] (last time)
587            # - also works in general when t == T[i]
588            ratio = 0
589        else:
590            #t is now between index and index+1
591            ratio = (t - self.time[self.index])/\
592                    (self.time[self.index+1] - self.time[self.index])
593
594        #Compute interpolated values
595        q = zeros(len(self.quantity_names), Float)
596        #print "self.precomputed_values", self.precomputed_values
597        for i, name in enumerate(self.quantity_names):
598            Q = self.precomputed_values[name]
599
600            if self.spatial is False:
601                #If there is no spatial info               
602                assert len(Q.shape) == 1
603
604                Q0 = Q[self.index]
605                if ratio > 0: Q1 = Q[self.index+1]
606
607            else:
608                if x is not None and y is not None:
609                    #Interpolate to x, y
610                   
611                    raise 'x,y interpolation not yet implemented'
612                else:
613                    #Use precomputed point
614                    Q0 = Q[self.index, point_id]
615                    if ratio > 0:
616                        Q1 = Q[self.index+1, point_id]
617           
618            #Linear temporal interpolation   
619            if ratio > 0:
620                if Q0 == NAN and Q1 == NAN:
621                    q[i]  = Q0
622                else:
623                    q[i] = Q0 + ratio*(Q1 - Q0)
624            else:
625                q[i] = Q0
626
627
628        #Return vector of interpolated values
629        #if len(q) == 1:
630        #    return q[0]
631        #else:
632        #    return q
633
634
635        #Return vector of interpolated values
636        #FIXME:
637        if self.spatial is True:
638            return q
639        else:
640            #Replicate q according to x and y
641            #This is e.g used for Wind_stress
642            if x is None or y is None: 
643                return q
644            else:
645                try:
646                    N = len(x)
647                except:
648                    return q
649                else:
650                    from Numeric import ones, Float
651                    #x is a vector - Create one constant column for each value
652                    N = len(x)
653                    assert len(y) == N, 'x and y must have same length'
654                    res = []
655                    for col in q:
656                        res.append(col*ones(N, Float))
657                       
658                return res
659
660
661    def get_time(self):
662        """Return model time as a vector of timesteps
663        """
664        return self.time
665
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
713
714def interpolate_sww(sww_file, time, interpolation_points,
715                    quantity_names=None, verbose=False):
716    """
717    obsolete.
718    use file_function in utils
719    """
720    #open sww file
721    x, y, volumes, time, quantities = read_sww(sww_file)
722    print "x",x
723    print "y",y
724   
725    print "time", time
726    print "quantities", quantities
727
728    #Add the x and y together
729    vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
730
731    #Will return the quantity values at the specified times and locations
732    interp = Interpolation_interface(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
804#-------------------------------------------------------------
805if __name__ == "__main__":
806    names = ["x","y"]
807    someiterable = [[1,2],[3,4]]
808    csvwriter = writer(file("some.csv", "wb"))
809    csvwriter.writerow(names)
810    for row in someiterable:
811        csvwriter.writerow(row)
Note: See TracBrowser for help on using the repository browser.