Python scipy.interpolate 模块，NearestNDInterpolator() 实例源码

def interp_plume_dv(plume_model,pts_file,fname_out='plume.xy',zero_center=True):
   '''
   Takes a plume model file, interpolates it to a new grid, and writes a new file

   args--------------------------------------------------------------------------
   plume_model: the plume model file path.  cartesian coordinate system
   pts_file: a file containing the points to interpolate to
   fname_out: the name of the output file
   '''

   f1 = np.loadtxt(plume_model) 
   x = f1[:,0]
   y = 6371.0 - f1[:,1] 
   v = f1[:,2]

   f2 = np.loadtxt(pts_file,skiprows=1)
   x_new = f2[:,0]
   y_new = f2[:,1]


   x_mirror = x*-1.0 
   x_pts = np.hstack((x,x_mirror))
   y_pts = np.hstack((y,y))
   v_pts = np.hstack((v,v))

   if not zero_center: 
      x_pts += max(x_new)/2

   interpolator = NearestNDInterpolator((x_pts,y_pts),v_pts)

   v_new = interpolator(x_new,y_new)

   f3 = open(fname_out,'w')
   for i in range(0,len(x_new)):
      f3.write('{} {} {}'.format(x_new[i],y_new[i],v_new[i])+'\n') 

   f3.close()

def _build_vertex_values(self, items, area_func, value_func):
        """
        Interpolate vertice with known altitudes to get altitudes for the remaining ones.
        """
        vertex_values = np.empty(self.vertices.shape[:1], dtype=np.int32)
        vertex_value_mask = np.full(self.vertices.shape[:1], fill_value=False, dtype=np.bool)

        for item in items:
            i_vertices = np.unique(self.faces[np.array(tuple(chain(*area_func(item).faces)))].flatten())
            vertex_values[i_vertices] = value_func(item, i_vertices)
            vertex_value_mask[i_vertices] = True

        if not np.all(vertex_value_mask):
            interpolate = NearestNDInterpolator(self.vertices[vertex_value_mask],
                                                vertex_values[vertex_value_mask])
            vertex_values[np.logical_not(vertex_value_mask)] = interpolate(
                *np.transpose(self.vertices[np.logical_not(vertex_value_mask)])
            )

        return vertex_values

def __init__(self, scenario: model.Scenario):
        """Initialize the model."""
        super().__init__(scenario)

        s = self._scenario

        flag_rs = flag_ss = flag_ns = 0
        if s.mechanism == 'SS':
            flag_ss = 1
        elif s.mechanism == 'NS':
            flag_ns = 1
        elif s.mechanism == 'RS':
            flag_rs = 1

        event = (s.mag, s.depth_hyp, flag_rs, flag_ss, flag_ns, s.dist_jb,
                 s.v_s30)

        global INTERPOLATOR
        if INTERPOLATOR is None:
            with np.load(fname_data) as data:
                INTERPOLATOR = NearestNDInterpolator(data['events'],
                                                     data['predictions'])
        prediction = INTERPOLATOR(event)
        self._ln_resp = prediction[0::2]
        self._ln_std = np.sqrt(prediction[1::2])

def surface_transform(vals, subj, hemi, direction="reverse"):
    """Transform a surface scalar map using spherical transform.

    Parameters
    ----------
    vals : array or Series
        Scalar value map to transform.
    subj : string
        Freesurfer subject ID.
    hemi : lh | rh
        Hemisphere data are defined on
    direction : reverse | forward
        Whether transformation should be from group space to subject space
        (reverse) or the other direction (forward).

    Returns
    -------
    out_vals : array or Series
        Scalar value map defined on new surface.

    """
    data_dir = PROJECT["data_dir"]
    sphere_reg_fname = op.join(data_dir, subj, "surf", hemi + ".sphere.reg")
    avg_sphere_fname = op.join(data_dir, "fsaverage/surf", hemi + ".sphere")

    sphere_reg, _ = nib.freesurfer.read_geometry(sphere_reg_fname)
    avg_sphere, _ = nib.freesurfer.read_geometry(avg_sphere_fname)

    if direction.startswith("f"):
        src_sphere, trg_sphere = sphere_reg, avg_sphere
    elif direction.startswith("r"):
        src_sphere, trg_sphere = avg_sphere, sphere_reg

    interpolator = interpolate.NearestNDInterpolator(src_sphere, vals)
    out_vals = interpolator(trg_sphere)

    if isinstance(vals, pd.Series):
        out_vals = pd.Series(out_vals)

    return out_vals

def fill_missing(mask, x, y, u, v):
    """Fill missing value with by interpolating the values from nearest
    neighbours"""
    # Construct the interpolators from the boundary values surrounding the
    # missing values.
    boundaries = find_boundaries(u.mask, mode='outer')
    points = np.stack((y[boundaries], x[boundaries])).T
    ip_u = NearestNDInterpolator(points, u[boundaries], rescale=False)
    ip_v = NearestNDInterpolator(points, v[boundaries], rescale=False)

    # Interpolate only missing values.
    missing = (y[mask], x[mask])
    u[mask] = ip_u(missing)
    v[mask] = ip_v(missing)

def __init__(self, data, categorical_parameters, continuous_parameters, func=None, order=1):
        self.key_columns = categorical_parameters
        self.parameter_columns = continuous_parameters
        self.func = func

        if len(self.parameter_columns) not in [1, 2]:
            raise ValueError("Only interpolation over 1 or 2 variables is supported")
        if len(self.parameter_columns) == 1 and order == 0:
            raise ValueError("Order 0 only supported for 2d interpolation")

        # These are the columns which the interpolation function will approximate
        value_columns = sorted(data.columns.difference(set(self.key_columns)|set(self.parameter_columns)))

        if self.key_columns:
            # Since there are key_columns we need to group the table by those
            # columns to get the sub-tables to fit
            sub_tables = data.groupby(self.key_columns)
        else:
            # There are no key columns so we will fit the whole table
            sub_tables = {None: data}.items()

        self.interpolations = {}

        for key, base_table in sub_tables:
            if base_table.empty:
                continue
            # For each permutation of the key columns build interpolations
            self.interpolations[key] = {}
            for value_column in value_columns:
                # For each value in the table build an interpolation function
                if len(self.parameter_columns) == 2:
                    # 2 variable interpolation
                    if order == 0:
                        x = base_table[list(self.parameter_columns)]
                        y = base_table[value_column]
                        func = interpolate.NearestNDInterpolator(x=x.values, y=y.values)
                    else:
                        index, column = self.parameter_columns
                        table = base_table.pivot(index=index, columns=column, values=value_column)
                        x = table.index.values
                        y = table.columns.values
                        z = table.values
                        func = interpolate.RectBivariateSpline(x=x, y=y, z=z, ky=order, kx=order).ev
                else:
                    # 1 variable interpolation
                    base_table = base_table.sort_values(by=self.parameter_columns[0])
                    x = base_table[self.parameter_columns[0]]
                    y = base_table[value_column]
                    func = interpolate.InterpolatedUnivariateSpline(x, y, k=order)
                self.interpolations[key][value_column] = func