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

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

Comments

File size: 42.6 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=vertex_coordinates,
87                                triangles=triangles,
88                                mesh_origin=mesh_origin,
89                                verbose=verbose,
90                                max_vertices_per_cell=max_vertices_per_cell)
91
92    def interpolate_polyline(self,
93                             f,
94                             vertex_coordinates,
95                             gauge_neighbour_id,
96                             point_coordinates=None,
97                             verbose=False):
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       
121        if isinstance(point_coordinates, Geospatial_data):
122            point_coordinates = point_coordinates.get_data_points( \
123                absolute = True)
124 
125        from utilities.polygon import point_on_line
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'
132        assert f.shape[0]==vertex_coordinates.shape[0], msg
133        msg='Must define function value at one or more nodes'
134        assert f.shape[0]>0, msg
135
136        n=f.shape[0]
137        if n==1:
138            z=f*z
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       
145        elif n>1:
146            for i in range(len(point_coordinates)):
147                found = False
148                for j in range(n):
149                    if gauge_neighbour_id[j]>=0:
150                        if point_on_line(point_coordinates[i],
151                                         [vertex_coordinates[j], vertex_coordinates[gauge_neighbour_id[j]]],
152                                         rtol=1.0e-6):
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                                 
167                if not found:
168                    z[i]=0.0
169                    #print 'point not on urs boundary'
170        return z
171
172    # FIXME: What is a good start_blocking_len value?
173    def interpolate(self,
174                    f,
175                    point_coordinates=None,
176                    start_blocking_len=500000,
177                    verbose=False):
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        """
201
202        # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the
203        # method is called even if interpolation points are unchanged.
204
205        #print "point_coordinates interpolate.interpolate", point_coordinates
206        if verbose: print 'Build intepolation object' 
207        if isinstance(point_coordinates, Geospatial_data):
208            point_coordinates = point_coordinates.get_data_points( \
209                absolute = True)
210
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 \
214                   len(self._point_coordinates) < start_blocking_len:
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)
226           
227        if point_coordinates is not None:
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:
235                #print 'BLOCKING'
236                #Handle blocking
237                self._A_can_be_reused = False
238                start = 0
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                   
248                for end in range(start_blocking_len,
249                                 len(point_coordinates),
250                                 start_blocking_len):
251                   
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
258                   
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
264   
265
266    def interpolate_block(self, f, point_coordinates, verbose=False):
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):
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
303        return self._get_point_data_z(f)
304   
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        """
313
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: 
320            z[i] = NAN
321        return z
322
323    def _build_interpolation_matrix_A(self,
324                                      point_coordinates,
325                                      verbose=False):
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
342        if verbose: print 'Building interpolation matrix'
343
344        # Convert point_coordinates to Numeric arrays, in case it was a list.
345        point_coordinates = ensure_numeric(point_coordinates, Float)
346       
347       
348        if verbose: print 'Getting indices inside mesh boundary'
349        self.inside_poly_indices, self.outside_poly_indices  = \
350                     in_and_outside_polygon(point_coordinates,
351                                            self.mesh.get_boundary_polygon(),
352                                            closed = True, verbose = verbose)
353       
354        #Build n x m interpolation matrix
355        if verbose and len(self.outside_poly_indices) > 0:
356            print '\n WARNING: Points outside mesh boundary. \n'
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           
361        m = self.mesh.number_of_nodes  # Nbr of basis functions (1/vertex)
362        n = point_coordinates.shape[0] # Nbr of data points
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
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)
375            x = point_coordinates[i]
376            element_found, sigma0, sigma1, sigma2, k = \
377                           search_tree_of_vertices(self.root, self.mesh, x)
378           
379            # Update interpolation matrix A if necessary
380            if element_found is True:
381                # Assign values to matrix A
382
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
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
397def benchmark_interpolate(vertices,
398                          vertex_attributes,
399                          triangles, points,
400                          max_points_per_cell=None,
401                          start_blocking_len=500000,
402                          mesh_origin=None):
403    """
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             
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, 
417                         max_vertices_per_cell=max_points_per_cell,
418                         mesh_origin=mesh_origin)
419           
420    calc = interp.interpolate(vertex_attributes
421                              ,points
422                              ,start_blocking_len=start_blocking_len)
423    #print "calc", calc
424   
425def interpolate_sww2csv(sww_file,
426                        points,
427                        depth_file,
428                        velocity_x_file,
429                        velocity_y_file,
430                        stage_file=None,
431                        #quantities = ['depth', 'velocity'],
432                        time_thinning=1,
433                        verbose=True,
434                        use_cache = True):
435    """
436    Interpolate the quantities at a given set of locations, given
437    an sww file.
438    The results are written to csv files.
439
440    sww_file is the input sww file.
441    points is a list of the 'gauges' x,y location.
442    depth_file is the name of the output depth file
443    velocity_x_file is the name of the output x velocity file.
444    velocity_y_file is the name of the output y velocity file.
445    stage_file is the name of the output stage file.
446
447    In the csv files columns represents the gauges and each row is a
448    time slice.
449   
450   
451    Time_thinning_number controls how many timesteps to use. Only
452        timesteps with index%time_thinning_number == 0 will used, or
453        in other words a value of 3, say, will cause the algorithm to
454        use every third time step.
455
456    In the future let points be a points file.
457    And let the user choose the quantities.
458
459    This is currently quite specific.
460    If it is need to be more general, change things.
461
462    """
463    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
464    #print "points",points
465    points = ensure_absolute(points)
466    point_count = len(points)
467    callable_sww = file_function(sww_file,
468                                 quantities=quantities,
469                                 interpolation_points=points,
470                                 verbose=verbose,
471                                 time_thinning=time_thinning,
472                                 use_cache=use_cache)
473   
474    depth_writer = writer(file(depth_file, "wb"))
475    velocity_x_writer = writer(file(velocity_x_file, "wb"))
476    velocity_y_writer = writer(file(velocity_y_file, "wb"))
477    if stage_file is not None:
478        stage_writer = writer(file(stage_file, "wb"))
479    # Write heading
480    heading = [str(x[0])+ ':' + str(x[1]) for x in points]
481    heading.insert(0, "time")
482    depth_writer.writerow(heading)
483    velocity_x_writer.writerow(heading)
484    velocity_y_writer.writerow(heading)
485    if stage_file is not None:
486        stage_writer.writerow(heading)   
487   
488    for time in callable_sww.get_time():
489        depths = [time]
490        velocity_xs = [time]
491        velocity_ys = [time]
492        if stage_file is not None: 
493            stages = [time] 
494        for point_i, point in enumerate(points):
495            quantities = callable_sww(time,point_i)
496            #print "quantities", quantities
497           
498            w = quantities[0]
499            z = quantities[1]
500            momentum_x = quantities[2]
501            momentum_y = quantities[3]
502            depth = w - z
503             
504            if w == NAN or z == NAN or momentum_x == NAN:
505                velocity_x = NAN
506            else:
507                if depth > 1.e-30: # use epsilon
508                    velocity_x = momentum_x / depth  #Absolute velocity
509                else:
510                    velocity_x = 0
511            if w == NAN or z == NAN or momentum_y == NAN:
512                velocity_y = NAN
513            else:
514                if depth > 1.e-30: # use epsilon
515                    velocity_y = momentum_y / depth  #Absolute velocity
516                else:
517                    velocity_y = 0
518            depths.append(depth)
519            velocity_xs.append(velocity_x)
520            velocity_ys.append(velocity_y)
521            if stage_file is not None:
522                stages.append(w)
523        depth_writer.writerow(depths)
524        velocity_x_writer.writerow(velocity_xs)
525        velocity_y_writer.writerow(velocity_ys)
526        if stage_file is not None:
527            stage_writer.writerow(stages)       
528
529
530class Interpolation_function:
531    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
532    which is interpolated from time series defined at vertices of
533    triangular mesh (such as those stored in sww files)
534
535    Let m be the number of vertices, n the number of triangles
536    and p the number of timesteps.
537    Also, let N be the number of interpolation points.
538
539    Mandatory input
540        time:               px1 array of monotonously increasing times (Float)
541        quantities:         Dictionary of arrays or 1 array (Float)
542                            The arrays must either have dimensions pxm or mx1.
543                            The resulting function will be time dependent in
544                            the former case while it will be constant with
545                            respect to time in the latter case.
546       
547    Optional input:
548        quantity_names:     List of keys into the quantities dictionary for
549                            imposing a particular order on the output vector.
550        vertex_coordinates: mx2 array of coordinates (Float)
551        triangles:          nx3 array of indices into vertex_coordinates (Int)
552        interpolation_points: Nx2 array of coordinates to be interpolated to
553        verbose:            Level of reporting
554   
555   
556    The quantities returned by the callable object are specified by
557    the list quantities which must contain the names of the
558    quantities to be returned and also reflect the order, e.g. for
559    the shallow water wave equation, on would have
560    quantities = ['stage', 'xmomentum', 'ymomentum']
561
562    The parameter interpolation_points decides at which points interpolated
563    quantities are to be computed whenever object is called.
564    If None, return average value
565
566    FIXME (Ole): Need to allow vertex coordinates and interpolation points to be
567    geospatial data objects
568
569    Time assumed to be relative to starttime (FIXME (Ole): This comment should be removed)
570    All coordinates assume origin of (0,0) - e.g. georeferencing must be taken care of
571    outside this function
572    """
573 
574   
575    def __init__(self,
576                 time,
577                 quantities,
578                 quantity_names=None, 
579                 vertex_coordinates=None,
580                 triangles=None,
581                 interpolation_points=None,
582                 time_thinning=1,
583                 verbose=False,
584                 gauge_neighbour_id=None):
585        """Initialise object and build spatial interpolation if required
586
587        Time_thinning_number controls how many timesteps to use. Only timesteps with
588        index%time_thinning_number == 0 will used, or in other words a value of 3, say,
589        will cause the algorithm to use every third time step.
590        """
591
592        from Numeric import array, zeros, Float, alltrue, concatenate,\
593             reshape, ArrayType
594
595
596        from anuga.config import time_format
597        import types
598
599
600        # Check temporal info
601        time = ensure_numeric(time)       
602        msg = 'Time must be a monotonuosly '
603        msg += 'increasing sequence %s' %time
604        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
605
606
607        # Check if quantities is a single array only
608        if type(quantities) != types.DictType:
609            quantities = ensure_numeric(quantities)
610            quantity_names = ['Attribute']
611
612            #Make it a dictionary
613            quantities = {quantity_names[0]: quantities}
614
615
616        # Use keys if no names are specified
617        if quantity_names is None:
618            quantity_names = quantities.keys()
619
620
621        # Check spatial info
622        if vertex_coordinates is None:
623            self.spatial = False
624        else:
625            # FIXME (Ole): Try ensure_numeric here -
626            #this function knows nothing about georefering.
627            vertex_coordinates = ensure_absolute(vertex_coordinates)
628
629            if triangles is not None:
630                triangles = ensure_numeric(triangles)
631            self.spatial = True         
632
633        # Thin timesteps if needed
634        # Note array() is used to make the thinned arrays contiguous in memory
635        self.time = array(time[::time_thinning])         
636        for name in quantity_names:
637            if len(quantities[name].shape) == 2:
638                quantities[name] = array(quantities[name][::time_thinning,:])
639             
640        # Save for use with statistics
641        self.quantities_range = {}
642        for name in quantity_names:
643            q = quantities[name][:].flat
644            self.quantities_range[name] = [min(q), max(q)]
645       
646        self.quantity_names = quantity_names       
647        self.vertex_coordinates = vertex_coordinates
648        self.interpolation_points = interpolation_points
649       
650
651        self.index = 0    # Initial time index
652        self.precomputed_values = {}
653       
654           
655        # Precomputed spatial interpolation if requested
656        if interpolation_points is not None:
657            #no longer true. sts files have spatial = True but
658            #if self.spatial is False:
659            #    raise 'Triangles and vertex_coordinates must be specified'
660            #
661            try:
662                self.interpolation_points = interpolation_points = ensure_numeric(interpolation_points)
663            except:
664                msg = 'Interpolation points must be an N x 2 Numeric array '+\
665                      'or a list of points\n'
666                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
667                                      '...')
668                raise msg
669
670            if triangles is not None and vertex_coordinates is not None:
671                # Check that all interpolation points fall within
672                # mesh boundary as defined by triangles and vertex_coordinates.
673                from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
674                from anuga.utilities.polygon import outside_polygon           
675
676                # Create temporary mesh object from mesh info passed
677                # into this function.
678                mesh = Mesh(vertex_coordinates, triangles)
679                mesh_boundary_polygon = mesh.get_boundary_polygon()
680
681           
682                indices = outside_polygon(interpolation_points,
683                                          mesh_boundary_polygon)
684
685                # Record result
686                #self.mesh_boundary_polygon = mesh_boundary_polygon
687                self.indices_outside_mesh = indices
688
689                # Report
690                if len(indices) > 0:
691                    msg = 'Interpolation points in Interpolation function fall ' 
692                    msg += 'outside specified mesh. '
693                    msg += 'Offending points:\n'
694                    out_interp_pts = []
695                    for i in indices:
696                        msg += '%d: %s\n' %(i, interpolation_points[i])
697                        out_interp_pts.append(ensure_numeric(interpolation_points[i]))
698
699                    if verbose is True:
700                        import sys
701                        if sys.platform == 'win32':
702                            from anuga.utilities.polygon import plot_polygons
703                            #out_interp_pts = take(interpolation_points,[indices])
704                            title = 'Interpolation points fall outside specified mesh'
705                            plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts],
706                                          ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose)
707
708                    # Joaquim Luis suggested this as an Exception, so
709                    # that the user can now what the problem is rather than
710                    # looking for NaN's. However, NANs are handy as they can
711                    # be ignored leaving good points for continued processing.
712                    if verbose:
713                        print msg
714                    #raise Exception(msg)
715            elif triangles is None and vertex_coordinates is not None:#jj
716                #Dealing with sts file
717                pass
718            else:
719                msg = 'Sww file function requires both triangles and vertex_coordinates. sts file file function requires the later.'
720                raise Exception(msg)
721
722            # Plot boundary and interpolation points
723            if verbose is True:
724                import sys
725                if sys.platform == 'win32':
726                    from anuga.utilities.polygon import plot_polygons
727                    title = 'Interpolation function: Polygon and interpolation points'
728                    plot_polygons([mesh_boundary_polygon,interpolation_points],
729                                         ['line','point'],figname='points_boundary',label=title,verbose=verbose)
730
731            m = len(self.interpolation_points)
732            p = len(self.time)
733           
734            for name in quantity_names:
735                self.precomputed_values[name] = zeros((p, m), Float)
736
737            # Build interpolator
738            if verbose:
739                if triangles is not None and vertex_coordinates is not None:
740                    msg = 'Building interpolation matrix from source mesh '
741                    msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
742                                                           triangles.shape[0])
743                elif triangles is None and vertex_coordinates is not None:
744                    msg = 'Building interpolation matrix from source points'
745               
746                print msg
747
748               
749            interpol = Interpolate(vertex_coordinates,
750                                   triangles,
751                                   verbose=verbose)
752
753            if verbose:
754                print 'Interpolating (%d interpolation points, %d timesteps).'\
755                      %(self.interpolation_points.shape[0], self.time.shape[0]),
756           
757                if time_thinning > 1:
758                    print 'Timesteps were thinned by a factor of %d' %time_thinning
759                else:
760                    print
761
762            for i, t in enumerate(self.time):
763                # Interpolate quantities at this timestep
764                if verbose and i%((p+10)/10)==0:
765                    print '  time step %d of %d' %(i, p)
766                   
767                for name in quantity_names:
768                    if len(quantities[name].shape) == 2:
769                        Q = quantities[name][i,:] # Quantities at timestep i
770                    else:
771                        Q = quantities[name][:]   # No time dependency
772
773                    if verbose and i%((p+10)/10)==0:
774                        print '    quantity %s, size=%d' %(name, len(Q))
775                       
776                    # Interpolate
777                    if triangles is not None and vertex_coordinates is not None:   
778                        result = interpol.interpolate(Q,
779                                                      point_coordinates=\
780                                                      self.interpolation_points,
781                                                      verbose=False) # Don't clutter
782                    elif triangles is None and vertex_coordinates is not None:
783                        result=interpol.interpolate_polyline(Q,vertex_coordinates,gauge_neighbour_id,point_coordinates=self.interpolation_points)
784
785                    #assert len(result), len(interpolation_points)
786                    self.precomputed_values[name][i, :] = result
787
788                   
789            # Report
790            if verbose:
791                print self.statistics()
792                #self.print_statistics()
793           
794        else:
795            # Store quantitites as is
796            for name in quantity_names:
797                self.precomputed_values[name] = quantities[name]
798
799    def __repr__(self):
800        # return 'Interpolation function (spatio-temporal)'
801        return self.statistics()
802 
803    def __call__(self, t, point_id=None, x=None, y=None):
804        """Evaluate f(t) or f(t, point_id)
805       
806        Inputs:
807          t: time - Model time. Must lie within existing timesteps
808          point_id: index of one of the preprocessed points.
809     
810                   
811          If spatial info is present and all of point_id
812          are None an exception is raised
813                   
814          If no spatial info is present, point_id arguments are ignored
815          making f a function of time only.
816
817
818          FIXME: f(t, x, y) x, y could overrided location, point_id ignored 
819          FIXME: point_id could also be a slice
820          FIXME: What if x and y are vectors?
821          FIXME: What about f(x,y) without t?
822        """
823
824        from math import pi, cos, sin, sqrt
825        from Numeric import zeros, Float
826        from anuga.abstract_2d_finite_volumes.util import mean       
827
828        if self.spatial is True:
829            if point_id is None:
830                if x is None or y is None:
831                    msg = 'Either point_id or x and y must be specified'
832                    raise Exception(msg)
833            else:
834                if self.interpolation_points is None:
835                    msg = 'Interpolation_function must be instantiated ' +\
836                          'with a list of interpolation points before parameter ' +\
837                          'point_id can be used'
838                    raise Exception(msg)
839
840        msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
841        msg += ' does not match model time: %.16f\n' %t
842        if t < self.time[0]: raise Exception(msg)
843        if t > self.time[-1]: raise Exception(msg)
844
845        oldindex = self.index #Time index
846
847        # Find current time slot
848        while t > self.time[self.index]: self.index += 1
849        while t < self.time[self.index]: self.index -= 1
850
851        if t == self.time[self.index]:
852            # Protect against case where t == T[-1] (last time)
853            #  - also works in general when t == T[i]
854            ratio = 0
855        else:
856            # t is now between index and index+1
857            ratio = (t - self.time[self.index])/\
858                    (self.time[self.index+1] - self.time[self.index])
859
860        # Compute interpolated values
861        q = zeros(len(self.quantity_names), Float)
862        # print "self.precomputed_values", self.precomputed_values
863        for i, name in enumerate(self.quantity_names):
864            Q = self.precomputed_values[name]
865
866            if self.spatial is False:
867                # If there is no spatial info               
868                assert len(Q.shape) == 1
869
870                Q0 = Q[self.index]
871                if ratio > 0: Q1 = Q[self.index+1]
872
873            else:
874                if x is not None and y is not None:
875                    # Interpolate to x, y
876                   
877                    raise 'x,y interpolation not yet implemented'
878                else:
879                    # Use precomputed point
880                    Q0 = Q[self.index, point_id]
881                    if ratio > 0:
882                        Q1 = Q[self.index+1, point_id]
883
884            # Linear temporal interpolation   
885            if ratio > 0:
886                if Q0 == NAN and Q1 == NAN:
887                    q[i] = Q0
888                else:
889                    q[i] = Q0 + ratio*(Q1 - Q0)
890            else:
891                q[i] = Q0
892
893
894        # Return vector of interpolated values
895        # if len(q) == 1:
896        #     return q[0]
897        # else:
898        #     return q
899
900
901        # Return vector of interpolated values
902        # FIXME:
903        if self.spatial is True:
904            return q
905        else:
906            # Replicate q according to x and y
907            # This is e.g used for Wind_stress
908            if x is None or y is None: 
909                return q
910            else:
911                try:
912                    N = len(x)
913                except:
914                    return q
915                else:
916                    from Numeric import ones, Float
917                    # x is a vector - Create one constant column for each value
918                    N = len(x)
919                    assert len(y) == N, 'x and y must have same length'
920                    res = []
921                    for col in q:
922                        res.append(col*ones(N, Float))
923                       
924                return res
925
926
927    def get_time(self):
928        """Return model time as a vector of timesteps
929        """
930        return self.time
931
932
933    def statistics(self):
934        """Output statistics about interpolation_function
935        """
936       
937        vertex_coordinates = self.vertex_coordinates
938        interpolation_points = self.interpolation_points               
939        quantity_names = self.quantity_names
940        #quantities = self.quantities
941        precomputed_values = self.precomputed_values                 
942               
943        x = vertex_coordinates[:,0]
944        y = vertex_coordinates[:,1]               
945
946        str =  '------------------------------------------------\n'
947        str += 'Interpolation_function (spatio-temporal) statistics:\n'
948        str += '  Extent:\n'
949        str += '    x in [%f, %f], len(x) == %d\n'\
950               %(min(x), max(x), len(x))
951        str += '    y in [%f, %f], len(y) == %d\n'\
952               %(min(y), max(y), len(y))
953        str += '    t in [%f, %f], len(t) == %d\n'\
954               %(min(self.time), max(self.time), len(self.time))
955        str += '  Quantities:\n'
956        for name in quantity_names:
957            minq, maxq = self.quantities_range[name]
958            str += '    %s in [%f, %f]\n' %(name, minq, maxq)           
959            #q = quantities[name][:].flat
960            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
961
962        if interpolation_points is not None:   
963            str += '  Interpolation points (xi, eta):'\
964                   ' number of points == %d\n' %interpolation_points.shape[0]
965            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
966                                            max(interpolation_points[:,0]))
967            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
968                                             max(interpolation_points[:,1]))
969            str += '  Interpolated quantities (over all timesteps):\n'
970       
971            for name in quantity_names:
972                q = precomputed_values[name][:].flat
973                str += '    %s at interpolation points in [%f, %f]\n'\
974                       %(name, min(q), max(q))
975        str += '------------------------------------------------\n'
976
977        return str
978
979
980def interpolate_sww(sww_file, time, interpolation_points,
981                    quantity_names=None, verbose=False):
982    """
983    obsolete.
984    use file_function in utils
985    """
986    #open sww file
987    x, y, volumes, time, quantities = read_sww(sww_file)
988    print "x",x
989    print "y",y
990   
991    print "time", time
992    print "quantities", quantities
993
994    #Add the x and y together
995    vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
996
997    #Will return the quantity values at the specified times and locations
998    interp = Interpolation_interface(time,
999                                     quantities,
1000                                     quantity_names=quantity_names, 
1001                                     vertex_coordinates=vertex_coordinates,
1002                                     triangles=volumes,
1003                                     interpolation_points=interpolation_points,
1004                                     verbose=verbose)
1005
1006
1007def read_sww(file_name):
1008    """
1009    obsolete - Nothing should be calling this
1010   
1011    Read in an sww file.
1012   
1013    Input;
1014    file_name - the sww file
1015   
1016    Output;
1017    x - Vector of x values
1018    y - Vector of y values
1019    z - Vector of bed elevation
1020    volumes - Array.  Each row has 3 values, representing
1021    the vertices that define the volume
1022    time - Vector of the times where there is stage information
1023    stage - array with respect to time and vertices (x,y)
1024    """
1025
1026    msg = 'Function read_sww in interpolat.py is obsolete'
1027    raise Exception, msg
1028   
1029    #FIXME Have this reader as part of data_manager?
1030   
1031    from Scientific.IO.NetCDF import NetCDFFile     
1032    import tempfile
1033    import sys
1034    import os
1035   
1036    #Check contents
1037    #Get NetCDF
1038   
1039    # see if the file is there.  Throw a QUIET IO error if it isn't
1040    # I don't think I could implement the above
1041   
1042    #throws prints to screen if file not present
1043    junk = tempfile.mktemp(".txt")
1044    fd = open(junk,'w')
1045    stdout = sys.stdout
1046    sys.stdout = fd
1047    fid = NetCDFFile(file_name, 'r') 
1048    sys.stdout = stdout
1049    fd.close()
1050    #clean up
1051    os.remove(junk)     
1052     
1053    # Get the variables
1054    x = fid.variables['x'][:]
1055    y = fid.variables['y'][:]
1056    volumes = fid.variables['volumes'][:] 
1057    time = fid.variables['time'][:]
1058
1059    keys = fid.variables.keys()
1060    keys.remove("x")
1061    keys.remove("y")
1062    keys.remove("volumes")
1063    keys.remove("time")
1064     #Turn NetCDF objects into Numeric arrays
1065    quantities = {}
1066    for name in keys:
1067        quantities[name] = fid.variables[name][:]
1068           
1069    fid.close()
1070    return x, y, volumes, time, quantities
1071
1072
1073#-------------------------------------------------------------
1074if __name__ == "__main__":
1075    names = ["x","y"]
1076    someiterable = [[1,2],[3,4]]
1077    csvwriter = writer(file("some.csv", "wb"))
1078    csvwriter.writerow(names)
1079    for row in someiterable:
1080        csvwriter.writerow(row)
Note: See TracBrowser for help on using the repository browser.