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

Last change on this file since 5566 was 5566, checked in by ole, 16 years ago

Manual updates from Ted Rigby.

File size: 43.2 KB
RevLine 
[2879]1"""Least squares interpolation.
2
3   Implements a least-squares interpolation.
[3409]4   Putting mesh data onto points.
[2879]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
[3319]13* The public interface to Interpolate
[2879]14__init__
15interpolate
16interpolate_block
17
18"""
19
20import time
21import os
22from warnings import warn
[3411]23from math import sqrt
24from csv import writer, DictWriter
[2879]25
26from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
27     ArrayType, allclose, take, NewAxis, arange
28
[3514]29from anuga.caching.caching import cache
[3560]30from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
[3514]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
[3768]36from anuga.geospatial_data.geospatial_data import Geospatial_data
37from anuga.geospatial_data.geospatial_data import ensure_absolute
[3514]38from anuga.fit_interpolate.search_functions import search_tree_of_vertices
39from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
[3560]40from anuga.abstract_2d_finite_volumes.util import file_function
[2879]41
42
43class Interpolate (FitInterpolate):
44       
45    def __init__(self,
46                 vertex_coordinates,
47                 triangles,
48                 mesh_origin=None,
49                 verbose=False,
[4614]50                 max_vertices_per_cell=None):
[2879]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        """
[3689]78
79        # FIXME (Ole): Need an input check
[2879]80       
81        # Initialise variabels
82        self._A_can_be_reused = False
83        self._point_coordinates = None
84       
85        FitInterpolate.__init__(self,
[4779]86                                vertex_coordinates=vertex_coordinates,
87                                triangles=triangles,
88                                mesh_origin=mesh_origin,
89                                verbose=verbose,
90                                max_vertices_per_cell=max_vertices_per_cell)
[2879]91
[5361]92    def interpolate_polyline(self,
93                             f,
94                             vertex_coordinates,
[5372]95                             gauge_neighbour_id,
[5361]96                             point_coordinates=None,
97                             verbose=False):
[5403]98        """Interpolate linearly between values f on nodes (vertex coordinates) of a polyline to midpoints of triangles
99        of boundary.
100
101        f is the data on the polyline nodes.
102
103        The mesh values representing a smooth surface are
104        assumed to be specified in f.
105
106        Inputs:
107          f: Vector or array of data at the polyline nodes.
108              If f is an array, interpolation will be done for each column as
109              per underlying matrix-matrix multiplication
110          point_coordinates: Interpolate polyline data to these positions.
111              List of coordinate pairs [x, y] of
112              data points or an nx2 Numeric array or a Geospatial_data object
113             
114        Output:
115          Interpolated values at inputted points (z).
116        """
117       
118
119        # FIXME: There is an option of passing a tolerance into this
120       
[5361]121        if isinstance(point_coordinates, Geospatial_data):
122            point_coordinates = point_coordinates.get_data_points( \
123                absolute = True)
124 
[5403]125        from utilities.polygon import point_on_line
[5361]126        from Numeric import ones
127        z=ones(len(point_coordinates),Float)
128
129        msg='point coordinates are not given (interpolate.py)'
130        assert point_coordinates is not None, msg
131        msg='function value must be specified at every interpolation node'
[5403]132        assert f.shape[0]==vertex_coordinates.shape[0], msg
[5361]133        msg='Must define function value at one or more nodes'
[5403]134        assert f.shape[0]>0, msg
[5361]135
136        n=f.shape[0]
137        if n==1:
138            z=f*z
[5403]139            msg = 'Polyline contained only one point. I need more. ', str(f)
140            raise Exception, msg
141           
142        # FIXME (John): add unit test for only One vertex point. Exception should be thrown.
143
144       
[5361]145        elif n>1:
146            for i in range(len(point_coordinates)):
147                found = False
[5372]148                for j in range(n):
149                    if gauge_neighbour_id[j]>=0:
[5403]150                        if point_on_line(point_coordinates[i],
151                                         [vertex_coordinates[j], vertex_coordinates[gauge_neighbour_id[j]]],
152                                         rtol=1.0e-6):
[5372]153                            found=True
154                            x0=vertex_coordinates[j][0]
155                            y0=vertex_coordinates[j][1]
156                            x1=vertex_coordinates[gauge_neighbour_id[j]][0]
157                            y1=vertex_coordinates[gauge_neighbour_id[j]][1]
158                            x2=point_coordinates[i][0]
159                            y2=point_coordinates[i][1]
160                           
161                            segment_len=sqrt((x1-x0)**2+(y1-y0)**2)
162                            dist=sqrt((x2-x0)**2+(y2-y0)**2)
163                            z[i]=(f[gauge_neighbour_id[j]]-f[j])/segment_len*dist+f[j]
164                            #print 'element found on segment'
165                            break
166                                 
[5361]167                if not found:
168                    z[i]=0.0
169                    #print 'point not on urs boundary'
170        return z
171
[3689]172    # FIXME: What is a good start_blocking_len value?
[3896]173    def interpolate(self,
174                    f,
175                    point_coordinates=None,
176                    start_blocking_len=500000,
177                    verbose=False):
[2879]178        """Interpolate mesh data f to determine values, z, at points.
179
180        f is the data on the mesh vertices.
181
182        The mesh values representing a smooth surface are
183        assumed to be specified in f.
184
185        Inputs:
186          f: Vector or array of data at the mesh vertices.
187              If f is an array, interpolation will be done for each column as
188              per underlying matrix-matrix multiplication
189          point_coordinates: Interpolate mesh data to these positions.
190              List of coordinate pairs [x, y] of
191              data points or an nx2 Numeric array or a Geospatial_data object
192             
193              If point_coordinates is absent, the points inputted last time
194              this method was called are used, if possible.
195          start_blocking_len: If the # of points is more or greater than this,
196              start blocking
197
198        Output:
199          Interpolated values at inputted points (z).
200        """
[3689]201
[3768]202        # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
203        # method is called even if interpolation points are unchanged.
[5222]204
[5321]205        #print "point_coordinates interpolate.interpolate", point_coordinates
206        if verbose: print 'Build intepolation object' 
[3689]207        if isinstance(point_coordinates, Geospatial_data):
[2879]208            point_coordinates = point_coordinates.get_data_points( \
209                absolute = True)
[3689]210
[2879]211        # Can I interpolate, based on previous point_coordinates?
212        if point_coordinates is None:
213            if self._A_can_be_reused is True and \
[3689]214                   len(self._point_coordinates) < start_blocking_len:
[2879]215                z = self._get_point_data_z(f,
216                                           verbose=verbose)
217            elif self._point_coordinates is not None:
218                #     if verbose, give warning
219                if verbose:
220                    print 'WARNING: Recalculating A matrix, due to blocking.'
221                point_coordinates = self._point_coordinates
222            else:
223                #There are no good point_coordinates. import sys; sys.exit()
224                msg = 'ERROR (interpolate.py): No point_coordinates inputted'
225                raise Exception(msg)
[5321]226           
227        if point_coordinates is not None:
[2879]228            self._point_coordinates = point_coordinates
229            if len(point_coordinates) < start_blocking_len or \
230                   start_blocking_len == 0:
231                self._A_can_be_reused = True
232                z = self.interpolate_block(f, point_coordinates,
233                                           verbose=verbose)
234            else:
[3689]235                #print 'BLOCKING'
[2879]236                #Handle blocking
237                self._A_can_be_reused = False
[3689]238                start = 0
[2879]239                # creating a dummy array to concatenate to.
240               
241                f = ensure_numeric(f, Float)
242                #print "f.shape",f.shape
243                if len(f.shape) > 1:
244                    z = zeros((0,f.shape[1]))
245                else:
246                    z = zeros((0,))
247                   
[3689]248                for end in range(start_blocking_len,
249                                 len(point_coordinates),
250                                 start_blocking_len):
[3768]251                   
[2879]252                    t = self.interpolate_block(f, point_coordinates[start:end],
253                                               verbose=verbose)
254                    #print "t", t
255                    #print "z", z
256                    z = concatenate((z,t))
257                    start = end
[3768]258                   
[2879]259                end = len(point_coordinates)
260                t = self.interpolate_block(f, point_coordinates[start:end],
261                                           verbose=verbose)
262                z = concatenate((z,t))
263        return z
[5222]264   
[2879]265
[5222]266    def interpolate_block(self, f, point_coordinates, verbose=False):
[2879]267        """
268        Call this if you want to control the blocking or make sure blocking
269        doesn't occur.
270
271        Return the point data, z.
272       
273        See interpolate for doc info.
274        """
275        if isinstance(point_coordinates,Geospatial_data):
[5222]276            point_coordinates = point_coordinates.get_data_points(\
277                absolute=True)
278
279        # Convert lists to Numeric arrays if necessary
280        point_coordinates = ensure_numeric(point_coordinates, Float)
281        f = ensure_numeric(f, Float)       
282           
283        self._A = self._build_interpolation_matrix_A(point_coordinates,
284                                                     verbose=verbose)
285
286
287        # Check that input dimensions are compatible
288        msg = 'Two colums must be specified in point coordinates. I got shape=%s'\
289              %(str(point_coordinates.shape))
290        assert point_coordinates.shape[1] == 2, msg
291
292        msg = 'The number of rows in matrix A must be the same as the number of points supplied.'
293        msg += ' I got %d points and %d matrix rows.'\
294               %(point_coordinates.shape[0], self._A.shape[0])
295        assert point_coordinates.shape[0] == self._A.shape[0], msg       
296
297        msg = 'The number of columns in matrix A must be the same as the number of mesh vertices.'
298        msg += ' I got %d vertices and %d matrix columns.'\
299               %(f.shape[0], self._A.shape[1])       
300        assert self._A.shape[1] == f.shape[0], msg
301
302        # Compute Matrix vector product and return
[2879]303        return self._get_point_data_z(f)
[5222]304   
[2879]305
306    def _get_point_data_z(self, f, verbose=False):
307        """
308        Return the point data, z.
309       
310        Precondition,
311        The _A matrix has been created
312        """
[5222]313
[2879]314        z = self._A * f
315        # Taking into account points outside the mesh.
316        #print "self.outside_poly_indices", self.outside_poly_indices
317        #print "self.inside_poly_indices", self.inside_poly_indices
318        #print "z", z
319        for i in self.outside_poly_indices: 
[3452]320            z[i] = NAN
[2879]321        return z
322
323    def _build_interpolation_matrix_A(self,
[3689]324                                      point_coordinates,
[3896]325                                      verbose=False):
[2879]326        """Build n x m interpolation matrix, where
327        n is the number of data points and
328        m is the number of basis functions phi_k (one per vertex)
329
330        This algorithm uses a quad tree data structure for fast binning
331        of data points
332        origin is a 3-tuple consisting of UTM zone, easting and northing.
333        If specified coordinates are assumed to be relative to this origin.
334
335        This one will override any data_origin that may be specified in
336        instance interpolation
337
338        Preconditions
339        Point_coordindates and mesh vertices have the same origin.
340        """
341
[5321]342        if verbose: print 'Building interpolation matrix'
[5222]343
344        # Convert point_coordinates to Numeric arrays, in case it was a list.
[2879]345        point_coordinates = ensure_numeric(point_coordinates, Float)
[4648]346       
[5222]347       
[2879]348        if verbose: print 'Getting indices inside mesh boundary'
349        self.inside_poly_indices, self.outside_poly_indices  = \
[3686]350                     in_and_outside_polygon(point_coordinates,
351                                            self.mesh.get_boundary_polygon(),
352                                            closed = True, verbose = verbose)
[4648]353       
[2879]354        #Build n x m interpolation matrix
[4589]355        if verbose and len(self.outside_poly_indices) > 0:
[3018]356            print '\n WARNING: Points outside mesh boundary. \n'
[3019]357        # Since you can block, throw a warning, not an error.
358        if verbose and 0 == len(self.inside_poly_indices):
359            print '\n WARNING: No points within the mesh! \n'
360           
[3945]361        m = self.mesh.number_of_nodes  # Nbr of basis functions (1/vertex)
362        n = point_coordinates.shape[0] # Nbr of data points
[2879]363
364        if verbose: print 'Number of datapoints: %d' %n
365        if verbose: print 'Number of basis functions: %d' %m
366
367        A = Sparse(n,m)
368
369        n = len(self.inside_poly_indices)
370        #Compute matrix elements for points inside the mesh
[4589]371        if verbose: print 'Building interpolation matrix from %d points' %n
372        for d, i in enumerate(self.inside_poly_indices):
373            # For each data_coordinate point
374            if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
[2879]375            x = point_coordinates[i]
376            element_found, sigma0, sigma1, sigma2, k = \
377                           search_tree_of_vertices(self.root, self.mesh, x)
[4589]378           
379            # Update interpolation matrix A if necessary
[2879]380            if element_found is True:
[4589]381                # Assign values to matrix A
[2879]382
[4589]383                j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0
384                j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
385                j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
[2879]386
387                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
388                js     = [j0,j1,j2]
389
390                for j in js:
391                    A[i,j] = sigmas[j]
392            else:
393                msg = 'Could not find triangle for point', x
394                raise Exception(msg)
395        return A
396
[4651]397def benchmark_interpolate(vertices,
398                          vertex_attributes,
399                          triangles, points,
400                          max_points_per_cell=None,
[4839]401                          start_blocking_len=500000,
402                          mesh_origin=None):
[4651]403    """
[4839]404    points: Interpolate mesh data to these positions.
405              List of coordinate pairs [x, y] of
406              data points or an nx2 Numeric array or a Geospatial_data object
407             
[4651]408    No test for this yet.
409    Note, this has no time the input data has no time dimension.  Which is
410    different from most of the data we interpolate, eg sww info.
411     
412        Output:
413          Interpolated values at inputted points.
414    """
415    interp = Interpolate(vertices,
416                         triangles, 
[4839]417                         max_vertices_per_cell=max_points_per_cell,
418                         mesh_origin=mesh_origin)
[4651]419           
420    calc = interp.interpolate(vertex_attributes
421                              ,points
422                              ,start_blocking_len=start_blocking_len)
[4839]423    #print "calc", calc
424   
[3419]425def interpolate_sww2csv(sww_file,
426                        points,
[3411]427                        depth_file,
[3456]428                        velocity_x_file,
429                        velocity_y_file,
[5196]430                        stage_file=None,
[5496]431                        froude_file=None,
[3417]432                        #quantities = ['depth', 'velocity'],
[5447]433                        time_thinning=1,
[3417]434                        verbose=True,
435                        use_cache = True):
[3411]436    """
437    Interpolate the quantities at a given set of locations, given
438    an sww file.
[5454]439    The results are written to csv files.
440
441    sww_file is the input sww file.
442    points is a list of the 'gauges' x,y location.
443    depth_file is the name of the output depth file
444    velocity_x_file is the name of the output x velocity file.
445    velocity_y_file is the name of the output y velocity file.
446    stage_file is the name of the output stage file.
447
448    In the csv files columns represents the gauges and each row is a
449    time slice.
[5447]450   
[5454]451   
[5447]452    Time_thinning_number controls how many timesteps to use. Only
453        timesteps with index%time_thinning_number == 0 will used, or
454        in other words a value of 3, say, will cause the algorithm to
455        use every third time step.
[3411]456
[4663]457    In the future let points be a points file.
[5454]458    And let the user choose the quantities.
[3417]459
460    This is currently quite specific.
[5454]461    If it is need to be more general, change things.
[3445]462
[3411]463    """
464    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
465    #print "points",points
466    points = ensure_absolute(points)
467    point_count = len(points)
468    callable_sww = file_function(sww_file,
469                                 quantities=quantities,
470                                 interpolation_points=points,
471                                 verbose=verbose,
[5447]472                                 time_thinning=time_thinning,
[3411]473                                 use_cache=use_cache)
474   
475    depth_writer = writer(file(depth_file, "wb"))
[3456]476    velocity_x_writer = writer(file(velocity_x_file, "wb"))
477    velocity_y_writer = writer(file(velocity_y_file, "wb"))
[5196]478    if stage_file is not None:
479        stage_writer = writer(file(stage_file, "wb"))
[5496]480    if froude_file is not None:
481        froude_writer = writer(file(froude_file, "wb"))
[3411]482    # Write heading
[3419]483    heading = [str(x[0])+ ':' + str(x[1]) for x in points]
484    heading.insert(0, "time")
485    depth_writer.writerow(heading)
[3456]486    velocity_x_writer.writerow(heading)
487    velocity_y_writer.writerow(heading)
[5196]488    if stage_file is not None:
[5496]489        stage_writer.writerow(heading) 
490    if froude_file is not None:
491        froude_writer.writerow(heading)     
[3456]492   
[3411]493    for time in callable_sww.get_time():
[3419]494        depths = [time]
[3456]495        velocity_xs = [time]
496        velocity_ys = [time]
[5196]497        if stage_file is not None: 
498            stages = [time] 
[5496]499        if froude_file is not None: 
500            froudes = [time] 
[3411]501        for point_i, point in enumerate(points):
502            quantities = callable_sww(time,point_i)
503            #print "quantities", quantities
504           
505            w = quantities[0]
506            z = quantities[1]
[3456]507            momentum_x = quantities[2]
508            momentum_y = quantities[3]
[3411]509            depth = w - z
510             
[3456]511            if w == NAN or z == NAN or momentum_x == NAN:
512                velocity_x = NAN
[3411]513            else:
514                if depth > 1.e-30: # use epsilon
[3456]515                    velocity_x = momentum_x / depth  #Absolute velocity
[3411]516                else:
[3456]517                    velocity_x = 0
518            if w == NAN or z == NAN or momentum_y == NAN:
519                velocity_y = NAN
520            else:
521                if depth > 1.e-30: # use epsilon
522                    velocity_y = momentum_y / depth  #Absolute velocity
523                else:
524                    velocity_y = 0
[5496]525            if depth < 1.e-30: # use epsilon
526                froude = NAN
527            else:
528                froude = sqrt(velocity_x*velocity_x + velocity_y*velocity_y)/ \
529                         sqrt(depth * 9.8066) # gravity m/s/s
[3411]530            depths.append(depth)
[3456]531            velocity_xs.append(velocity_x)
532            velocity_ys.append(velocity_y)
[5196]533            if stage_file is not None:
534                stages.append(w)
[5496]535            if froude_file is not None:
536                froudes.append(froude)
[3411]537        depth_writer.writerow(depths)
[3456]538        velocity_x_writer.writerow(velocity_xs)
539        velocity_y_writer.writerow(velocity_ys)
[5196]540        if stage_file is not None:
[5496]541            stage_writer.writerow(stages)   
542        if froude_file is not None:
543            froude_writer.writerow(froudes)           
[3411]544
[3768]545
[2879]546class Interpolation_function:
547    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
548    which is interpolated from time series defined at vertices of
549    triangular mesh (such as those stored in sww files)
550
551    Let m be the number of vertices, n the number of triangles
[5288]552    and p the number of timesteps.
553    Also, let N be the number of interpolation points.
[2879]554
555    Mandatory input
556        time:               px1 array of monotonously increasing times (Float)
557        quantities:         Dictionary of arrays or 1 array (Float)
558                            The arrays must either have dimensions pxm or mx1.
559                            The resulting function will be time dependent in
[5288]560                            the former case while it will be constant with
[2879]561                            respect to time in the latter case.
562       
563    Optional input:
[5288]564        quantity_names:     List of keys into the quantities dictionary for
565                            imposing a particular order on the output vector.
[2879]566        vertex_coordinates: mx2 array of coordinates (Float)
567        triangles:          nx3 array of indices into vertex_coordinates (Int)
568        interpolation_points: Nx2 array of coordinates to be interpolated to
569        verbose:            Level of reporting
570   
571   
572    The quantities returned by the callable object are specified by
573    the list quantities which must contain the names of the
574    quantities to be returned and also reflect the order, e.g. for
575    the shallow water wave equation, on would have
576    quantities = ['stage', 'xmomentum', 'ymomentum']
577
578    The parameter interpolation_points decides at which points interpolated
579    quantities are to be computed whenever object is called.
580    If None, return average value
[3896]581
582    FIXME (Ole): Need to allow vertex coordinates and interpolation points to be
583    geospatial data objects
[3900]584
[5288]585    Time assumed to be relative to starttime (FIXME (Ole): This comment should be removed)
586    All coordinates assume origin of (0,0) - e.g. georeferencing must be taken care of
587    outside this function
[2879]588    """
[3768]589 
[2879]590   
591    def __init__(self,
592                 time,
593                 quantities,
[3850]594                 quantity_names=None, 
595                 vertex_coordinates=None,
596                 triangles=None,
597                 interpolation_points=None,
[3900]598                 time_thinning=1,
[5372]599                 verbose=False,
600                 gauge_neighbour_id=None):
[2879]601        """Initialise object and build spatial interpolation if required
[3900]602
603        Time_thinning_number controls how many timesteps to use. Only timesteps with
604        index%time_thinning_number == 0 will used, or in other words a value of 3, say,
605        will cause the algorithm to use every third time step.
[2879]606        """
607
608        from Numeric import array, zeros, Float, alltrue, concatenate,\
609             reshape, ArrayType
610
611
[3514]612        from anuga.config import time_format
[2879]613        import types
614
[4396]615
[3850]616        # Check temporal info
[2879]617        time = ensure_numeric(time)       
618        msg = 'Time must be a monotonuosly '
619        msg += 'increasing sequence %s' %time
620        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
621
622
[3850]623        # Check if quantities is a single array only
[2879]624        if type(quantities) != types.DictType:
625            quantities = ensure_numeric(quantities)
626            quantity_names = ['Attribute']
627
[5566]628            # Make it a dictionary
[2879]629            quantities = {quantity_names[0]: quantities}
630
631
[3850]632        # Use keys if no names are specified
[2879]633        if quantity_names is None:
634            quantity_names = quantities.keys()
635
636
[3850]637        # Check spatial info
[2879]638        if vertex_coordinates is None:
639            self.spatial = False
[5288]640        else:
641            # FIXME (Ole): Try ensure_numeric here -
642            #this function knows nothing about georefering.
[3941]643            vertex_coordinates = ensure_absolute(vertex_coordinates)
[2879]644
[5361]645            if triangles is not None:
646                triangles = ensure_numeric(triangles)
647            self.spatial = True         
[2879]648
[3900]649        # Thin timesteps if needed
650        # Note array() is used to make the thinned arrays contiguous in memory
651        self.time = array(time[::time_thinning])         
652        for name in quantity_names:
653            if len(quantities[name].shape) == 2:
654                quantities[name] = array(quantities[name][::time_thinning,:])
[2879]655             
[3850]656        # Save for use with statistics
[3050]657        self.quantities_range = {}
658        for name in quantity_names:
659            q = quantities[name][:].flat
660            self.quantities_range[name] = [min(q), max(q)]
661       
[2879]662        self.quantity_names = quantity_names       
663        self.vertex_coordinates = vertex_coordinates
664        self.interpolation_points = interpolation_points
[3900]665       
666
[2879]667        self.index = 0    # Initial time index
668        self.precomputed_values = {}
[3050]669       
[4396]670           
[3850]671        # Precomputed spatial interpolation if requested
[2879]672        if interpolation_points is not None:
[5361]673            #no longer true. sts files have spatial = True but
674            #if self.spatial is False:
675            #    raise 'Triangles and vertex_coordinates must be specified'
676            #
[2879]677            try:
[4833]678                self.interpolation_points = interpolation_points = ensure_numeric(interpolation_points)
[2879]679            except:
680                msg = 'Interpolation points must be an N x 2 Numeric array '+\
681                      'or a list of points\n'
682                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
683                                      '...')
684                raise msg
685
[5361]686            if triangles is not None and vertex_coordinates is not None:
687                # Check that all interpolation points fall within
688                # mesh boundary as defined by triangles and vertex_coordinates.
689                from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
690                from anuga.utilities.polygon import outside_polygon           
[2879]691
[5361]692                # Create temporary mesh object from mesh info passed
693                # into this function.
694                mesh = Mesh(vertex_coordinates, triangles)
695                mesh_boundary_polygon = mesh.get_boundary_polygon()
[4833]696
[4820]697           
[5361]698                indices = outside_polygon(interpolation_points,
699                                          mesh_boundary_polygon)
[4820]700
[5361]701                # Record result
702                #self.mesh_boundary_polygon = mesh_boundary_polygon
703                self.indices_outside_mesh = indices
[4820]704
[5361]705                # Report
706                if len(indices) > 0:
707                    msg = 'Interpolation points in Interpolation function fall ' 
708                    msg += 'outside specified mesh. '
709                    msg += 'Offending points:\n'
710                    out_interp_pts = []
711                    for i in indices:
712                        msg += '%d: %s\n' %(i, interpolation_points[i])
713                        out_interp_pts.append(ensure_numeric(interpolation_points[i]))
[4820]714
[5361]715                    if verbose is True:
716                        import sys
717                        if sys.platform == 'win32':
718                            from anuga.utilities.polygon import plot_polygons
719                            #out_interp_pts = take(interpolation_points,[indices])
720                            title = 'Interpolation points fall outside specified mesh'
721                            plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts],
722                                          ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose)
[4851]723
[5361]724                    # Joaquim Luis suggested this as an Exception, so
725                    # that the user can now what the problem is rather than
726                    # looking for NaN's. However, NANs are handy as they can
727                    # be ignored leaving good points for continued processing.
728                    if verbose:
729                        print msg
730                    #raise Exception(msg)
731            elif triangles is None and vertex_coordinates is not None:#jj
732                #Dealing with sts file
733                pass
734            else:
735                msg = 'Sww file function requires both triangles and vertex_coordinates. sts file file function requires the later.'
736                raise Exception(msg)
[4820]737
[4833]738            # Plot boundary and interpolation points
[4845]739            if verbose is True:
[4869]740                import sys
741                if sys.platform == 'win32':
742                    from anuga.utilities.polygon import plot_polygons
743                    title = 'Interpolation function: Polygon and interpolation points'
744                    plot_polygons([mesh_boundary_polygon,interpolation_points],
745                                         ['line','point'],figname='points_boundary',label=title,verbose=verbose)
[4820]746
[2879]747            m = len(self.interpolation_points)
[2884]748            p = len(self.time)
[4396]749           
[2879]750            for name in quantity_names:
751                self.precomputed_values[name] = zeros((p, m), Float)
752
[3850]753            # Build interpolator
[3900]754            if verbose:
[5361]755                if triangles is not None and vertex_coordinates is not None:
756                    msg = 'Building interpolation matrix from source mesh '
757                    msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
758                                                           triangles.shape[0])
759                elif triangles is None and vertex_coordinates is not None:
760                    msg = 'Building interpolation matrix from source points'
761               
[3900]762                print msg
763
764               
[2879]765            interpol = Interpolate(vertex_coordinates,
[3689]766                                   triangles,
[3850]767                                   verbose=verbose)
[2879]768
[3900]769            if verbose:
770                print 'Interpolating (%d interpolation points, %d timesteps).'\
771                      %(self.interpolation_points.shape[0], self.time.shape[0]),
772           
773                if time_thinning > 1:
774                    print 'Timesteps were thinned by a factor of %d' %time_thinning
775                else:
776                    print
777
[2884]778            for i, t in enumerate(self.time):
[3768]779                # Interpolate quantities at this timestep
[2879]780                if verbose and i%((p+10)/10)==0:
[5222]781                    print '  time step %d of %d' %(i, p)
[2879]782                   
783                for name in quantity_names:
784                    if len(quantities[name].shape) == 2:
[3768]785                        Q = quantities[name][i,:] # Quantities at timestep i
[2879]786                    else:
[3768]787                        Q = quantities[name][:]   # No time dependency
[5222]788
789                    if verbose and i%((p+10)/10)==0:
790                        print '    quantity %s, size=%d' %(name, len(Q))
791                       
[5361]792                    # Interpolate
793                    if triangles is not None and vertex_coordinates is not None:   
794                        result = interpol.interpolate(Q,
795                                                      point_coordinates=\
796                                                      self.interpolation_points,
797                                                      verbose=False) # Don't clutter
798                    elif triangles is None and vertex_coordinates is not None:
[5372]799                        result=interpol.interpolate_polyline(Q,vertex_coordinates,gauge_neighbour_id,point_coordinates=self.interpolation_points)
[4588]800
801                    #assert len(result), len(interpolation_points)
[2879]802                    self.precomputed_values[name][i, :] = result
[3768]803
[2879]804                   
[3768]805            # Report
[2879]806            if verbose:
807                print self.statistics()
808                #self.print_statistics()
809           
810        else:
[3768]811            # Store quantitites as is
812            for name in quantity_names:
[2879]813                self.precomputed_values[name] = quantities[name]
814
815    def __repr__(self):
[3768]816        # return 'Interpolation function (spatio-temporal)'
[2879]817        return self.statistics()
[3768]818 
819    def __call__(self, t, point_id=None, x=None, y=None):
[4370]820        """Evaluate f(t) or f(t, point_id)
821       
[2879]822        Inputs:
823          t: time - Model time. Must lie within existing timesteps
824          point_id: index of one of the preprocessed points.
[4370]825     
826                   
827          If spatial info is present and all of point_id
[2879]828          are None an exception is raised
829                   
[4370]830          If no spatial info is present, point_id arguments are ignored
[2879]831          making f a function of time only.
832
[4370]833
[4548]834          FIXME: f(t, x, y) x, y could overrided location, point_id ignored 
[2879]835          FIXME: point_id could also be a slice
836          FIXME: What if x and y are vectors?
837          FIXME: What about f(x,y) without t?
838        """
839
840        from math import pi, cos, sin, sqrt
841        from Numeric import zeros, Float
[3560]842        from anuga.abstract_2d_finite_volumes.util import mean       
[2879]843
844        if self.spatial is True:
845            if point_id is None:
846                if x is None or y is None:
847                    msg = 'Either point_id or x and y must be specified'
848                    raise Exception(msg)
849            else:
850                if self.interpolation_points is None:
851                    msg = 'Interpolation_function must be instantiated ' +\
852                          'with a list of interpolation points before parameter ' +\
853                          'point_id can be used'
854                    raise Exception(msg)
855
[2884]856        msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
[2879]857        msg += ' does not match model time: %.16f\n' %t
[2884]858        if t < self.time[0]: raise Exception(msg)
859        if t > self.time[-1]: raise Exception(msg)
[2879]860
861        oldindex = self.index #Time index
862
[4820]863        # Find current time slot
[2884]864        while t > self.time[self.index]: self.index += 1
865        while t < self.time[self.index]: self.index -= 1
[2879]866
[2884]867        if t == self.time[self.index]:
[4820]868            # Protect against case where t == T[-1] (last time)
869            #  - also works in general when t == T[i]
[2879]870            ratio = 0
871        else:
[4820]872            # t is now between index and index+1
[2884]873            ratio = (t - self.time[self.index])/\
874                    (self.time[self.index+1] - self.time[self.index])
[2879]875
[4820]876        # Compute interpolated values
[2879]877        q = zeros(len(self.quantity_names), Float)
[4820]878        # print "self.precomputed_values", self.precomputed_values
[2879]879        for i, name in enumerate(self.quantity_names):
880            Q = self.precomputed_values[name]
881
882            if self.spatial is False:
[4820]883                # If there is no spatial info               
[2879]884                assert len(Q.shape) == 1
885
886                Q0 = Q[self.index]
887                if ratio > 0: Q1 = Q[self.index+1]
888
889            else:
890                if x is not None and y is not None:
[4820]891                    # Interpolate to x, y
[2879]892                   
893                    raise 'x,y interpolation not yet implemented'
894                else:
[4820]895                    # Use precomputed point
[2879]896                    Q0 = Q[self.index, point_id]
897                    if ratio > 0:
898                        Q1 = Q[self.index+1, point_id]
[4588]899
[4820]900            # Linear temporal interpolation   
[2879]901            if ratio > 0:
[3452]902                if Q0 == NAN and Q1 == NAN:
[5482]903                    q[i] = Q0
[2879]904                else:
905                    q[i] = Q0 + ratio*(Q1 - Q0)
906            else:
907                q[i] = Q0
908
909
[4820]910        # Return vector of interpolated values
911        # if len(q) == 1:
912        #     return q[0]
913        # else:
914        #     return q
[2879]915
916
[4820]917        # Return vector of interpolated values
918        # FIXME:
[2879]919        if self.spatial is True:
920            return q
921        else:
[4820]922            # Replicate q according to x and y
923            # This is e.g used for Wind_stress
[2879]924            if x is None or y is None: 
925                return q
926            else:
927                try:
928                    N = len(x)
929                except:
930                    return q
931                else:
932                    from Numeric import ones, Float
[4820]933                    # x is a vector - Create one constant column for each value
[2879]934                    N = len(x)
935                    assert len(y) == N, 'x and y must have same length'
936                    res = []
937                    for col in q:
938                        res.append(col*ones(N, Float))
939                       
940                return res
941
942
[2884]943    def get_time(self):
944        """Return model time as a vector of timesteps
945        """
946        return self.time
947
[3768]948
[2879]949    def statistics(self):
950        """Output statistics about interpolation_function
951        """
952       
953        vertex_coordinates = self.vertex_coordinates
954        interpolation_points = self.interpolation_points               
955        quantity_names = self.quantity_names
[3055]956        #quantities = self.quantities
[2879]957        precomputed_values = self.precomputed_values                 
958               
959        x = vertex_coordinates[:,0]
960        y = vertex_coordinates[:,1]               
961
962        str =  '------------------------------------------------\n'
963        str += 'Interpolation_function (spatio-temporal) statistics:\n'
964        str += '  Extent:\n'
965        str += '    x in [%f, %f], len(x) == %d\n'\
966               %(min(x), max(x), len(x))
967        str += '    y in [%f, %f], len(y) == %d\n'\
968               %(min(y), max(y), len(y))
969        str += '    t in [%f, %f], len(t) == %d\n'\
[2884]970               %(min(self.time), max(self.time), len(self.time))
[2879]971        str += '  Quantities:\n'
972        for name in quantity_names:
[3050]973            minq, maxq = self.quantities_range[name]
974            str += '    %s in [%f, %f]\n' %(name, minq, maxq)           
975            #q = quantities[name][:].flat
976            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
[2879]977
978        if interpolation_points is not None:   
979            str += '  Interpolation points (xi, eta):'\
980                   ' number of points == %d\n' %interpolation_points.shape[0]
981            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
982                                            max(interpolation_points[:,0]))
983            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
984                                             max(interpolation_points[:,1]))
985            str += '  Interpolated quantities (over all timesteps):\n'
986       
987            for name in quantity_names:
988                q = precomputed_values[name][:].flat
989                str += '    %s at interpolation points in [%f, %f]\n'\
990                       %(name, min(q), max(q))
991        str += '------------------------------------------------\n'
992
993        return str
994
[3768]995
[2879]996def interpolate_sww(sww_file, time, interpolation_points,
[3896]997                    quantity_names=None, verbose=False):
[2879]998    """
999    obsolete.
1000    use file_function in utils
1001    """
1002    #open sww file
1003    x, y, volumes, time, quantities = read_sww(sww_file)
1004    print "x",x
1005    print "y",y
1006   
1007    print "time", time
1008    print "quantities", quantities
1009
1010    #Add the x and y together
1011    vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
1012
1013    #Will return the quantity values at the specified times and locations
[3768]1014    interp = Interpolation_interface(time,
1015                                     quantities,
1016                                     quantity_names=quantity_names, 
1017                                     vertex_coordinates=vertex_coordinates,
1018                                     triangles=volumes,
1019                                     interpolation_points=interpolation_points,
1020                                     verbose=verbose)
[2879]1021
1022
1023def read_sww(file_name):
1024    """
1025    obsolete - Nothing should be calling this
1026   
1027    Read in an sww file.
1028   
1029    Input;
1030    file_name - the sww file
1031   
1032    Output;
1033    x - Vector of x values
1034    y - Vector of y values
1035    z - Vector of bed elevation
1036    volumes - Array.  Each row has 3 values, representing
1037    the vertices that define the volume
1038    time - Vector of the times where there is stage information
1039    stage - array with respect to time and vertices (x,y)
1040    """
[4548]1041
1042    msg = 'Function read_sww in interpolat.py is obsolete'
1043    raise Exception, msg
[2879]1044   
1045    #FIXME Have this reader as part of data_manager?
1046   
1047    from Scientific.IO.NetCDF import NetCDFFile     
1048    import tempfile
1049    import sys
1050    import os
1051   
1052    #Check contents
1053    #Get NetCDF
1054   
1055    # see if the file is there.  Throw a QUIET IO error if it isn't
1056    # I don't think I could implement the above
1057   
1058    #throws prints to screen if file not present
1059    junk = tempfile.mktemp(".txt")
1060    fd = open(junk,'w')
1061    stdout = sys.stdout
1062    sys.stdout = fd
1063    fid = NetCDFFile(file_name, 'r') 
1064    sys.stdout = stdout
1065    fd.close()
1066    #clean up
1067    os.remove(junk)     
1068     
1069    # Get the variables
1070    x = fid.variables['x'][:]
1071    y = fid.variables['y'][:]
1072    volumes = fid.variables['volumes'][:] 
1073    time = fid.variables['time'][:]
1074
1075    keys = fid.variables.keys()
1076    keys.remove("x")
1077    keys.remove("y")
1078    keys.remove("volumes")
1079    keys.remove("time")
1080     #Turn NetCDF objects into Numeric arrays
1081    quantities = {}
1082    for name in keys:
1083        quantities[name] = fid.variables[name][:]
1084           
1085    fid.close()
1086    return x, y, volumes, time, quantities
1087
[3768]1088
[2879]1089#-------------------------------------------------------------
1090if __name__ == "__main__":
[3411]1091    names = ["x","y"]
1092    someiterable = [[1,2],[3,4]]
1093    csvwriter = writer(file("some.csv", "wb"))
1094    csvwriter.writerow(names)
1095    for row in someiterable:
1096        csvwriter.writerow(row)
Note: See TracBrowser for help on using the repository browser.