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

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

Introduced a covariance measure to verify Okushiri timeseries plus
some incidental work

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