Changeset 7673


Ignore:
Timestamp:
Mar 31, 2010, 10:57:09 PM (15 years ago)
Author:
James Hudson
Message:

Ticket 113 - added an output_centroids parameter to allow centroid data to be written to gauges, rather than the data at the sampled point.

Added 2 new unit tests for this functionality.

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py

    r7672 r7673  
    180180                                     verbose=verbose,
    181181                                     use_cache=use_cache,
    182                                                                         output_centroids = output_centroids)
     182                                    output_centroids = output_centroids)
    183183
    184184        if quake_offset_time is None:
    … …  
    190190                quake_time = time + quake_offset_time
    191191                points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug
    192                 point_quantities = callable_sww(time,point_i)
    193                
     192                point_quantities = callable_sww(time,point_i) # __call__ is overridden
     193                               
    194194                for quantity in quantities:
    195195                    if quantity == NAN:
    … …  
    271271                   title_on=None,
    272272                   use_cache=False,
    273                    verbose=False):
     273                   verbose=False,
     274                                   output_centroids=False):
    274275    """ Read sww file and plot the time series for the
    275276    prescribed quantities at defined gauge locations and
    … …  
    388389                        title_on,
    389390                        use_cache,
    390                         verbose)
     391                        verbose,
     392                                                output_centroids = output_centroids)
    391393    return k
    392394
    … …  
    423425                    title_on = None,
    424426                    use_cache = False,
    425                     verbose = False):   
     427                    verbose = False,
     428                                        output_centroids = False):   
    426429       
    427430    # FIXME(Ole): Shouldn't print statements here be governed by verbose?
    … …  
    486489                          time_thinning = time_thinning,
    487490                          verbose = verbose,
    488                           use_cache = use_cache)
     491                          use_cache = use_cache,
     492                                                  output_centroids = output_centroids)
    489493
    490494        # determine which gauges are contained in sww file

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py

    r7672 r7673  
    451451
    452452        # create a csv file containing our gauge points
    453         points = [[2.0,1.0],[5.0,5.0]]
     453        points = [[2.0,1.0],[4.5,4.0]]
    454454
    455455        points_file = tempfile.mktemp(".csv")
    456456        file_id = open(points_file,"w")
    457 # These values are where the centroids should be
     457# These values are where the centroids should be               
    458458#        file_id.write("name, easting, northing, elevation \n\
    459 #point1, 1.5, 1.5, 3.0\n\
    460 #point2, 4.5, 4.5, 9.0\n")
     459#point1, 2.0, 2.0, 3.0\n\
     460#point2, 4.0, 4.0, 9.0\n")
    461461 
    462  # These values are slightly off the centroids - will it find them?
     462# These values are slightly off the centroids - will it find the centroids?
    463463        file_id.write("name, easting, northing, elevation \n\
    464464point1, 2.0, 1.0, 3.0\n\
    465 point2, 5.5, 4.0, 9.0\n")
     465point2, 4.5, 4.0, 9.0\n")
    466466
    467467 
    … …  
    474474                       output_centroids=True)
    475475
    476         point1_answers_array = [[0.0,0.0,1.0,2.5,-1.5,3.0,4.0], [2.0,2.0/3600.,10.0,11.5,-1.5,3.0,4.0]]
     476        point1_answers_array = [[0.0,0.0,1.0,3.0,-2.0,3.0,4.0], [2.0,2.0/3600.,10.0,12.0,-2.0,3.0,4.0]]
    477477        point1_filename = 'gauge_point1.csv'
    478478        point1_handle = file(point1_filename)
    … …  
    482482        line=[]
    483483        for i,row in enumerate(point1_reader):
    484 #           print 'i',i,'row',row
    485484            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
    486485                         float(row[4]),float(row[5]),float(row[6])])
    … …  
    488487            assert num.allclose(line[i], point1_answers_array[i])
    489488
    490         point2_answers_array = [[0.0,0.0,1.0,5.5,-4.5,3.0,4.0], [2.0,2.0/3600.,10.0,14.5,-4.5,3.0,4.0]]
     489        point2_answers_array = [[0.0,0.0,1.0,5.0,-4.0,3.0,4.0], [2.0,2.0/3600.,10.0,14.0,-4.0,3.0,4.0]]
    491490        point2_filename = 'gauge_point2.csv'
    492491        point2_handle = file(point2_filename)
    … …  
    496495        line=[]
    497496        for i,row in enumerate(point2_reader):
    498 #           print 'i',i,'row',row
    499497            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
    500498                         float(row[4]),float(row[5]),float(row[6])])
    501 #           print 'assert line',line[i],'point1',point1_answers_array[i]
     499#           print i, 'assert line',line[i],'point2',point2_answers_array[i]
    502500            assert num.allclose(line[i], point2_answers_array[i])
    503501                         

  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r7672 r7673  
    5454                  use_cache=False,
    5555                  boundary_polygon=None,
    56                                   output_centroids=False):
     56                  output_centroids=False):
    5757    """Read time history of spatial data from NetCDF file and return
    5858    a callable object.
    … …  
    140140              'verbose': verbose,
    141141              'boundary_polygon': boundary_polygon,
    142                           'output_centroids': output_centroids}
     142              'output_centroids': output_centroids}
    143143
    144144    # Call underlying engine with or without caching
    … …  
    232232                                        verbose=verbose,
    233233                                        boundary_polygon=boundary_polygon,
    234                                                                                 output_centroids=output_centroids)
     234                                        output_centroids=output_centroids)
    235235    else:
    236236        # FIXME (Ole): Could add csv file here to address Ted Rigby's
    … …  
    260260                             verbose=False,
    261261                             boundary_polygon=None,
    262                                                         output_centroids=False):
     262                            output_centroids=False):
    263263    """Read time history of spatial data from NetCDF sww file and
    264264    return a callable object f(t,x,y)
    … …  
    501501                                   time_thinning=time_thinning,
    502502                                   verbose=verbose,
    503                                    gauge_neighbour_id=gauge_neighbour_id),
     503                                   gauge_neighbour_id=gauge_neighbour_id,
     504                                                                   output_centroids=output_centroids),
    504505            starttime)
    505506

  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r7326 r7673  
    6969                start_blocking_len=500000,
    7070                use_cache=False,
    71                 verbose=False):
     71                verbose=False,
     72                                output_centroids=False):
    7273    """Interpolate vertex_values to interpolation points.
    7374
    … …  
    153154    result = I.interpolate_block(vertex_values, interpolation_points,
    154155                                 use_cache=use_cache,
    155                                  verbose=verbose)
     156                                 verbose=verbose,
     157                                                                 output_centroids=output_centroids)
    156158
    157159    return result
    … …  
    228230                    point_coordinates=None,
    229231                    start_blocking_len=500000,
    230                     verbose=False):
     232                    verbose=False,
     233                                        output_centroids=False):
    231234        """Interpolate mesh data f to determine values, z, at points.
    232235
    … …  
    288291                self._A_can_be_reused = True
    289292                z = self.interpolate_block(f, point_coordinates,
    290                                            verbose=verbose)
     293                                           verbose=verbose, output_centroids=output_centroids)
    291294            else:
    292295                # Handle blocking
    … …  
    305308                                 start_blocking_len):
    306309                    t = self.interpolate_block(f, point_coordinates[start:end],
    307                                                verbose=verbose)
     310                                               verbose=verbose, output_centroids=output_centroids)
    308311                    z = num.concatenate((z, t), axis=0)    #??default#
    309312                    start = end
    … …  
    311314                end = len(point_coordinates)
    312315                t = self.interpolate_block(f, point_coordinates[start:end],
    313                                            verbose=verbose)
     316                                           verbose=verbose, output_centroids=output_centroids)
    314317                z = num.concatenate((z, t), axis=0)    #??default#
    315318        return z
    … …  
    317320
    318321    ##
    319     # @brief Control whether blocking occurs or not.
    320     # @param f ??
    321     # @param point_coordinates ??
    322     # @param use_cache ??
     322    # @brief Interpolate a block of vertices
     323    # @param f Array of arbitrary data to be interpolated
     324    # @param point_coordinates List of vertices to intersect with the mesh
     325    # @param use_cache True if caching should be used to accelerate the calculations
    323326    # @param verbose True if this function is verbose.
    324327    # @return ??
    325328    def interpolate_block(self, f, point_coordinates,
    326                           use_cache=False, verbose=False):
     329                          use_cache=False, verbose=False, output_centroids=False):
    327330        """
    328331        Call this if you want to control the blocking or make sure blocking
    … …  
    354357
    355358                X = cache(self._build_interpolation_matrix_A,
    356                           args=(point_coordinates,),
     359                          args=(point_coordinates, output_centroids),
    357360                          kwargs={'verbose': verbose},
    358361                          verbose=verbose)
    … …  
    372375                if reuse_A is False:
    373376                    X = self._build_interpolation_matrix_A(point_coordinates,
     377                                                           output_centroids,
    374378                                                           verbose=verbose)
    375379                    self.interpolation_matrices[key] = (X, point_coordinates)
    376380        else:
    377             X = self._build_interpolation_matrix_A(point_coordinates,
     381            X = self._build_interpolation_matrix_A(point_coordinates, output_centroids,
    378382                                                   verbose=verbose)
    379383
    … …  
    403407
    404408    ##
    405     # @brief ??
    406     # @param f ??
     409    # @brief Get interpolated data at given points.
     410        #        Applies a transform to all points to calculate the
     411        #        interpolated values. Points outside the mesh are returned as NaN.
     412        # @note self._A matrix must be valid
     413    # @param f Array of arbitrary data
    407414    # @param verbose True if this function is to be verbose.
    408     # @return ??
     415    # @return f transformed by interpolation matrix (f')
    409416    def _get_point_data_z(self, f, verbose=False):
    410417        """
    … …  
    424431    ##
    425432    # @brief Build NxM interpolation matrix.
    426     # @param point_coordinates ??
     433    # @param point_coordinates Points to sample at
     434        # @param output_centroids set to True to always sample from the centre
     435        #                         of the intersected triangle, instead of the intersection
     436        #                         point.
    427437    # @param verbose True if this function is to be verbose.
    428     # @return ??
     438    # @return Interpolation matrix A, plus lists of the points inside and outside the mesh.
    429439    def _build_interpolation_matrix_A(self,
    430440                                      point_coordinates,
     441                                      output_centroids=False,
    431442                                      verbose=False):
    432443        """Build n x m interpolation matrix, where
    … …  
    495506                j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1
    496507                j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
    497 
    498                 sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
    499                 js     = [j0, j1, j2]
    500 
    501                 for j in js:
    502                     A[i, j] = sigmas[j]
     508                js = [j0, j1, j2]
     509                               
     510                if output_centroids is False:
     511                    # Weight each vertex according to its distance from x
     512                    sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
     513                    for j in js:
     514                        A[i, j] = sigmas[j]
     515                else:
     516                                    # If centroids are needed, weight all 3 vertices equally
     517                    for j in js:
     518                        A[i, j] = 1.0/3.0                   
    503519            else:
    504520                msg = 'Could not find triangle for point', x
    … …  
    751767                 time_thinning=1,
    752768                 verbose=False,
    753                  gauge_neighbour_id=None):
     769                 gauge_neighbour_id=None,
     770                                 output_centroids=False):
    754771        """Initialise object and build spatial interpolation if required
    755772
    … …  
    976993                                                      point_coordinates=\
    977994                                                      self.interpolation_points,
    978                                                       verbose=False) # No clutter
     995                                                      verbose=False,
     996                                                                                                          output_centroids=output_centroids)
    979997                    elif triangles is None and vertex_coordinates is not None:
    980998                        result = interpolate_polyline(Q,
    … …  
    10021020
    10031021    ##
    1004     # @brief Override object() method.
     1022    # @brief Evaluate interpolation function
    10051023    # @param t Model time - must lie within existing timesteps.
    10061024    # @param point_id Index of one of the preprocessed points.

  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r7342 r7673  
    391391        assert num.allclose(num.sum(results, axis=1), 1.0)
    392392
     393               
     394    def test_arbitrary_datapoints_return_centroids(self):
     395        #Try arbitrary datapoints, making sure they return
     396        #an interpolation matrix for the intersected triangle's
     397        #centroid.
     398       
     399        a = [1.0, 0.0]
     400        b = [0.0, 3.0]
     401        c = [4.0,0.0]
     402        points = [a, b, c]
     403        vertices = [ [1,0,2] ]
     404
     405        data = [ [1.2, 1.5], [1.123, 1.768], [2.43, 0.44] ]
     406
     407        interp = Interpolate(points, vertices)
     408       
     409        third = [1.0/3.0, 1.0/3.0, 1.0/3.0]
     410        answer = [third, third, third]
     411               
     412        A,_,_ = interp._build_interpolation_matrix_A(data, output_centroids=True)
     413        results = A.todense()
     414        assert num.allclose(results, answer)           
     415               
     416               
    393417    def test_arbitrary_datapoints_some_outside(self):
    394418        #Try arbitrary datapoints one outside the triangle.
Note: See TracChangeset for help on using the changeset viewer.