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

Last change on this file since 5447 was 5447, checked in by duncan, 16 years ago

Current Hinwood scenario

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