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

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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