source: branches/inundation-numpy-branch/fit_interpolate/interpolate.py @ 7248

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