source: inundation/fit_interpolate/interpolate.py @ 3452

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

In numerical_tools, changing INF name to NAN.

File size: 29.2 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 caching.caching import cache
30from pyvolution.neighbour_mesh import Mesh
31from utilities.sparse import Sparse, Sparse_CSR
32from utilities.cg_solve import conjugate_gradient, VectorShapeError
33from coordinate_transforms.geo_reference import Geo_reference
34from pyvolution.quad import build_quadtree
35from utilities.numerical_tools import ensure_numeric, mean, NAN
36from utilities.polygon import in_and_outside_polygon
37from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute
38from fit_interpolate.search_functions import search_tree_of_vertices
39from fit_interpolate.general_fit_interpolate import FitInterpolate
40from 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_file,
284                        #quantities = ['depth', 'velocity'],
285                        verbose=True,
286                        use_cache = True):
287    """
288    Interpolate the quantities at a given set of locations, given
289    an sww file.
290    The results are written to a csv file.
291
292    In the future let points be a csv or xya file.
293    And the user choose the quantities.
294
295    This is currently quite specific.
296    If it need to be more general, chagne things.
297
298    This is really returning speed, not velocity.
299    """
300    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
301    #print "points",points
302    points = ensure_absolute(points)
303    point_count = len(points)
304    callable_sww = file_function(sww_file,
305                                 quantities=quantities,
306                                 interpolation_points=points,
307                                 verbose=verbose,
308                                 use_cache=use_cache)
309   
310    depth_writer = writer(file(depth_file, "wb"))
311    velocity_writer = writer(file(velocity_file, "wb"))
312    # Write heading
313    heading = [str(x[0])+ ':' + str(x[1]) for x in points]
314    heading.insert(0, "time")
315    depth_writer.writerow(heading)
316    velocity_writer.writerow(heading)
317    #for row in someiterable:
318        #csvwriter.writerow(row)
319    for time in callable_sww.get_time():
320        depths = [time]
321        velocitys = [time]
322        for point_i, point in enumerate(points):
323            quantities = callable_sww(time,point_i)
324            #print "quantities", quantities
325           
326            w = quantities[0]
327            z = quantities[1]
328            uh = quantities[2]
329            vh = quantities[3]
330            depth = w - z
331             
332            if w == NAN or z == NAN or uh == NAN or vh == NAN:
333                velocity = NAN
334            else:
335                momentum = sqrt(uh*uh + vh*vh)
336                if depth > 1.e-30: # use epsilon
337                    velocity = momentum / depth  #Absolute velocity
338                else:
339                    velocity = 0
340            depths.append(depth)
341            velocitys.append(velocity)
342        depth_writer.writerow(depths)
343        velocity_writer.writerow(velocitys)
344
345class Interpolation_function:
346    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
347    which is interpolated from time series defined at vertices of
348    triangular mesh (such as those stored in sww files)
349
350    Let m be the number of vertices, n the number of triangles
351    and p the number of timesteps.
352
353    Mandatory input
354        time:               px1 array of monotonously increasing times (Float)
355        quantities:         Dictionary of arrays or 1 array (Float)
356                            The arrays must either have dimensions pxm or mx1.
357                            The resulting function will be time dependent in
358                            the former case while it will be constan with
359                            respect to time in the latter case.
360       
361    Optional input:
362        quantity_names:     List of keys into the quantities dictionary
363        vertex_coordinates: mx2 array of coordinates (Float)
364        triangles:          nx3 array of indices into vertex_coordinates (Int)
365        interpolation_points: Nx2 array of coordinates to be interpolated to
366        verbose:            Level of reporting
367   
368   
369    The quantities returned by the callable object are specified by
370    the list quantities which must contain the names of the
371    quantities to be returned and also reflect the order, e.g. for
372    the shallow water wave equation, on would have
373    quantities = ['stage', 'xmomentum', 'ymomentum']
374
375    The parameter interpolation_points decides at which points interpolated
376    quantities are to be computed whenever object is called.
377    If None, return average value
378    """
379
380   
381   
382    def __init__(self,
383                 time,
384                 quantities,
385                 quantity_names = None, 
386                 vertex_coordinates = None,
387                 triangles = None,
388                 interpolation_points = None,
389                 verbose = False):
390        """Initialise object and build spatial interpolation if required
391        """
392
393        from Numeric import array, zeros, Float, alltrue, concatenate,\
394             reshape, ArrayType
395
396
397        #from util import mean, ensure_numeric
398        from config import time_format
399        import types
400
401
402        #Check temporal info
403        time = ensure_numeric(time)       
404        msg = 'Time must be a monotonuosly '
405        msg += 'increasing sequence %s' %time
406        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
407
408
409        #Check if quantities is a single array only
410        if type(quantities) != types.DictType:
411            quantities = ensure_numeric(quantities)
412            quantity_names = ['Attribute']
413
414            #Make it a dictionary
415            quantities = {quantity_names[0]: quantities}
416
417
418        #Use keys if no names are specified
419        if quantity_names is None:
420            quantity_names = quantities.keys()
421
422
423        #Check spatial info
424        if vertex_coordinates is None:
425            self.spatial = False
426        else:   
427            vertex_coordinates = ensure_numeric(vertex_coordinates)
428
429            assert triangles is not None, 'Triangles array must be specified'
430            triangles = ensure_numeric(triangles)
431            self.spatial = True           
432
433             
434        #Save for use with statistics
435
436        self.quantities_range = {}
437        for name in quantity_names:
438            q = quantities[name][:].flat
439            self.quantities_range[name] = [min(q), max(q)]
440       
441        self.quantity_names = quantity_names       
442        #self.quantities = quantities  #Takes too much memory     
443        self.vertex_coordinates = vertex_coordinates
444        self.interpolation_points = interpolation_points
445        self.time = time[:]  # Time assumed to be relative to starttime
446        self.index = 0    # Initial time index
447        self.precomputed_values = {}
448       
449
450
451
452       
453           
454        #Precomputed spatial interpolation if requested
455        if interpolation_points is not None:
456            if self.spatial is False:
457                raise 'Triangles and vertex_coordinates must be specified'
458           
459            try:
460                self.interpolation_points = ensure_numeric(interpolation_points)
461            except:
462                msg = 'Interpolation points must be an N x 2 Numeric array '+\
463                      'or a list of points\n'
464                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
465                                      '...')
466                raise msg
467
468
469            m = len(self.interpolation_points)
470            p = len(self.time)
471           
472            for name in quantity_names:
473                self.precomputed_values[name] = zeros((p, m), Float)
474
475            #Build interpolator
476            interpol = Interpolate(vertex_coordinates,
477                                     triangles,
478                                     #point_coordinates = \
479                                     #self.interpolation_points,
480                                     #alpha = 0,
481                                     verbose = verbose)
482
483            if verbose: print 'Interpolate'
484            for i, t in enumerate(self.time):
485                #Interpolate quantities at this timestep
486                if verbose and i%((p+10)/10)==0:
487                    print ' time step %d of %d' %(i, p)
488                   
489                for name in quantity_names:
490                    if len(quantities[name].shape) == 2:
491                        result = interpol.interpolate(quantities[name][i,:],
492                                     point_coordinates = \
493                                     self.interpolation_points)
494                    else:
495                       #Assume no time dependency
496                       result = interpol.interpolate(quantities[name][:],
497                                     point_coordinates = \
498                                     self.interpolation_points)
499                       
500                    self.precomputed_values[name][i, :] = result
501                   
502            #Report
503            if verbose:
504                print self.statistics()
505                #self.print_statistics()
506           
507        else:
508            #Store quantitites as is
509            for name in quantity_names:
510                self.precomputed_values[name] = quantities[name]
511
512
513    def __repr__(self):
514        #return 'Interpolation function (spatio-temporal)'
515        return self.statistics()
516   
517
518    def __call__(self, t, point_id = None, x = None, y = None):
519        """Evaluate f(t), f(t, point_id) or f(t, x, y)
520
521        Inputs:
522          t: time - Model time. Must lie within existing timesteps
523          point_id: index of one of the preprocessed points.
524          x, y:     Overrides location, point_id ignored
525         
526          If spatial info is present and all of x,y,point_id
527          are None an exception is raised
528                   
529          If no spatial info is present, point_id and x,y arguments are ignored
530          making f a function of time only.
531
532         
533          FIXME: point_id could also be a slice
534          FIXME: What if x and y are vectors?
535          FIXME: What about f(x,y) without t?
536        """
537
538        from math import pi, cos, sin, sqrt
539        from Numeric import zeros, Float
540        from util import mean       
541
542        if self.spatial is True:
543            if point_id is None:
544                if x is None or y is None:
545                    msg = 'Either point_id or x and y must be specified'
546                    raise Exception(msg)
547            else:
548                if self.interpolation_points is None:
549                    msg = 'Interpolation_function must be instantiated ' +\
550                          'with a list of interpolation points before parameter ' +\
551                          'point_id can be used'
552                    raise Exception(msg)
553
554        msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
555        msg += ' does not match model time: %.16f\n' %t
556        if t < self.time[0]: raise Exception(msg)
557        if t > self.time[-1]: raise Exception(msg)
558
559        oldindex = self.index #Time index
560
561        #Find current time slot
562        while t > self.time[self.index]: self.index += 1
563        while t < self.time[self.index]: self.index -= 1
564
565        if t == self.time[self.index]:
566            #Protect against case where t == T[-1] (last time)
567            # - also works in general when t == T[i]
568            ratio = 0
569        else:
570            #t is now between index and index+1
571            ratio = (t - self.time[self.index])/\
572                    (self.time[self.index+1] - self.time[self.index])
573
574        #Compute interpolated values
575        q = zeros(len(self.quantity_names), Float)
576        #print "self.precomputed_values", self.precomputed_values
577        for i, name in enumerate(self.quantity_names):
578            Q = self.precomputed_values[name]
579
580            if self.spatial is False:
581                #If there is no spatial info               
582                assert len(Q.shape) == 1
583
584                Q0 = Q[self.index]
585                if ratio > 0: Q1 = Q[self.index+1]
586
587            else:
588                if x is not None and y is not None:
589                    #Interpolate to x, y
590                   
591                    raise 'x,y interpolation not yet implemented'
592                else:
593                    #Use precomputed point
594                    Q0 = Q[self.index, point_id]
595                    if ratio > 0:
596                        Q1 = Q[self.index+1, point_id]
597           
598            #Linear temporal interpolation   
599            if ratio > 0:
600                if Q0 == NAN and Q1 == NAN:
601                    q[i]  = Q0
602                else:
603                    q[i] = Q0 + ratio*(Q1 - Q0)
604            else:
605                q[i] = Q0
606
607
608        #Return vector of interpolated values
609        #if len(q) == 1:
610        #    return q[0]
611        #else:
612        #    return q
613
614
615        #Return vector of interpolated values
616        #FIXME:
617        if self.spatial is True:
618            return q
619        else:
620            #Replicate q according to x and y
621            #This is e.g used for Wind_stress
622            if x is None or y is None: 
623                return q
624            else:
625                try:
626                    N = len(x)
627                except:
628                    return q
629                else:
630                    from Numeric import ones, Float
631                    #x is a vector - Create one constant column for each value
632                    N = len(x)
633                    assert len(y) == N, 'x and y must have same length'
634                    res = []
635                    for col in q:
636                        res.append(col*ones(N, Float))
637                       
638                return res
639
640
641    def get_time(self):
642        """Return model time as a vector of timesteps
643        """
644        return self.time
645
646    def statistics(self):
647        """Output statistics about interpolation_function
648        """
649       
650        vertex_coordinates = self.vertex_coordinates
651        interpolation_points = self.interpolation_points               
652        quantity_names = self.quantity_names
653        #quantities = self.quantities
654        precomputed_values = self.precomputed_values                 
655               
656        x = vertex_coordinates[:,0]
657        y = vertex_coordinates[:,1]               
658
659        str =  '------------------------------------------------\n'
660        str += 'Interpolation_function (spatio-temporal) statistics:\n'
661        str += '  Extent:\n'
662        str += '    x in [%f, %f], len(x) == %d\n'\
663               %(min(x), max(x), len(x))
664        str += '    y in [%f, %f], len(y) == %d\n'\
665               %(min(y), max(y), len(y))
666        str += '    t in [%f, %f], len(t) == %d\n'\
667               %(min(self.time), max(self.time), len(self.time))
668        str += '  Quantities:\n'
669        for name in quantity_names:
670            minq, maxq = self.quantities_range[name]
671            str += '    %s in [%f, %f]\n' %(name, minq, maxq)           
672            #q = quantities[name][:].flat
673            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
674
675        if interpolation_points is not None:   
676            str += '  Interpolation points (xi, eta):'\
677                   ' number of points == %d\n' %interpolation_points.shape[0]
678            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
679                                            max(interpolation_points[:,0]))
680            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
681                                             max(interpolation_points[:,1]))
682            str += '  Interpolated quantities (over all timesteps):\n'
683       
684            for name in quantity_names:
685                q = precomputed_values[name][:].flat
686                str += '    %s at interpolation points in [%f, %f]\n'\
687                       %(name, min(q), max(q))
688        str += '------------------------------------------------\n'
689
690        return str
691
692def interpolate_sww(sww_file, time, interpolation_points,
693                    quantity_names = None, verbose = False):
694    """
695    obsolete.
696    use file_function in utils
697    """
698    #open sww file
699    x, y, volumes, time, quantities = read_sww(sww_file)
700    print "x",x
701    print "y",y
702   
703    print "time", time
704    print "quantities", quantities
705
706    #Add the x and y together
707    vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
708
709    #Will return the quantity values at the specified times and locations
710    interp =  Interpolation_interface(
711                 time,
712                 quantities,
713                 quantity_names = quantity_names, 
714                 vertex_coordinates = vertex_coordinates,
715                 triangles = volumes,
716                 interpolation_points = interpolation_points,
717                 verbose = verbose)
718
719
720def read_sww(file_name):
721    """
722    obsolete - Nothing should be calling this
723   
724    Read in an sww file.
725   
726    Input;
727    file_name - the sww file
728   
729    Output;
730    x - Vector of x values
731    y - Vector of y values
732    z - Vector of bed elevation
733    volumes - Array.  Each row has 3 values, representing
734    the vertices that define the volume
735    time - Vector of the times where there is stage information
736    stage - array with respect to time and vertices (x,y)
737    """
738   
739    #FIXME Have this reader as part of data_manager?
740   
741    from Scientific.IO.NetCDF import NetCDFFile     
742    import tempfile
743    import sys
744    import os
745   
746    #Check contents
747    #Get NetCDF
748   
749    # see if the file is there.  Throw a QUIET IO error if it isn't
750    # I don't think I could implement the above
751   
752    #throws prints to screen if file not present
753    junk = tempfile.mktemp(".txt")
754    fd = open(junk,'w')
755    stdout = sys.stdout
756    sys.stdout = fd
757    fid = NetCDFFile(file_name, 'r') 
758    sys.stdout = stdout
759    fd.close()
760    #clean up
761    os.remove(junk)     
762     
763    # Get the variables
764    x = fid.variables['x'][:]
765    y = fid.variables['y'][:]
766    volumes = fid.variables['volumes'][:] 
767    time = fid.variables['time'][:]
768
769    keys = fid.variables.keys()
770    keys.remove("x")
771    keys.remove("y")
772    keys.remove("volumes")
773    keys.remove("time")
774     #Turn NetCDF objects into Numeric arrays
775    quantities = {}
776    for name in keys:
777        quantities[name] = fid.variables[name][:]
778           
779    fid.close()
780    return x, y, volumes, time, quantities
781
782#-------------------------------------------------------------
783if __name__ == "__main__":
784    names = ["x","y"]
785    someiterable = [[1,2],[3,4]]
786    csvwriter = writer(file("some.csv", "wb"))
787    csvwriter.writerow(names)
788    for row in someiterable:
789        csvwriter.writerow(row)
Note: See TracBrowser for help on using the repository browser.