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

Last change on this file since 4648 was 4648, checked in by duncan, 17 years ago

comments and updating benchmark

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