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

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

Refactored get_flow_through_cross_section in preparation for
get_energy_through_cross_section. These function work only on sww files.
Runtime versions still need to be done.

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