Matlab转Python：非默认坐标系下网格剖面速度提取问题

Got it, let's break down how to replicate Matlab's improfile function in Python using NumPy and SciPy, tailored to your custom grid setup where East/North are in the same coordinate system as your X/Y grid.

What Matlab's improfile does here #

Your Matlab line [cE_U, cN_U, c_U] = improfile(X,Y,U,East,North,100); does two key things:

1. Generates 100 evenly spaced sample points along the polyline defined by (East[0], North[0]) → (East[1], North[1]) → ...
2. Interpolates the velocity data U (on the X/Y grid) at each of those sample points, returning the sample coordinates and interpolated values.

Step 1: Generate evenly spaced points along the East/North polyline #

First, we need to create the 100 sample points along your custom polyline. We'll calculate the cumulative distance along the line, then interpolate to get uniform samples:

import numpy as np
from scipy.interpolate import interp1d

# Convert East and North to numpy arrays (if not already)
East = np.array(East)
North = np.array(North)

# Calculate cumulative distance along the polyline
dx = np.diff(East)
dy = np.diff(North)
segment_lengths = np.sqrt(dx**2 + dy**2)
cumulative_length = np.concatenate([[0], np.cumsum(segment_lengths)])
total_length = cumulative_length[-1]

# Generate 100 evenly spaced distance points
sample_distances = np.linspace(0, total_length, 100)

# Interpolate to get sample coordinates (cE_U, cN_U)
interp_east = interp1d(cumulative_length, East, kind='linear')
interp_north = interp1d(cumulative_length, North, kind='linear')
cE_U = interp_east(sample_distances)
cN_U = interp_north(sample_distances)

Step 2: Interpolate velocity U at the sample points #

Next, we'll interpolate the 40×40 velocity data U at our sample points. Since your X/Y define a custom grid, we'll use either RegularGridInterpolator (for regular grids, faster) or griddata (works for any grid):

Option 1: Using RegularGridInterpolator (best for regular grids)

If X and Y are 1D arrays defining the grid (e.g., from np.meshgrid), or can be easily extracted to 1D:

from scipy.interpolate import RegularGridInterpolator

# If X/Y are 2D (from meshgrid), extract the 1D axes:
# X_1d = X[0, :]
# Y_1d = Y[:, 0]

# Create the interpolator (note: U's shape should match (len(Y_1d), len(X_1d)))
interp_vel = RegularGridInterpolator((Y_1d, X_1d), U, method='linear')

# Interpolate at sample points (cE_U, cN_U)
# We need to pass points as (north, east) to match the grid axes order
sample_points = np.column_stack((cN_U, cE_U))
c_U = interp_vel(sample_points)

Option 2: Using griddata (works for any grid, regular or irregular)

If your X/Y grid is irregular, use this approach:

from scipy.interpolate import griddata

# Flatten the grid coordinates and velocity values
grid_points = np.column_stack((X.ravel(), Y.ravel()))
vel_values = U.ravel()

# Interpolate at sample points
c_U = griddata(grid_points, vel_values, (cE_U, cN_U), method='linear')

Putting it all together #

Combine the steps into a function that mirrors Matlab's improfile behavior:

def improfile(X, Y, U, East, North, n_samples):
    # Step 1: Generate sample points along polyline
    East = np.array(East)
    North = np.array(North)
    dx = np.diff(East)
    dy = np.diff(North)
    segment_lengths = np.sqrt(dx**2 + dy**2)
    cumulative_length = np.concatenate([[0], np.cumsum(segment_lengths)])
    sample_distances = np.linspace(0, cumulative_length[-1], n_samples)
    interp_east = interp1d(cumulative_length, East, kind='linear')
    interp_north = interp1d(cumulative_length, North, kind='linear')
    cE_U = interp_east(sample_distances)
    cN_U = interp_north(sample_distances)
    
    # Step 2: Interpolate velocity
    if X.ndim == 1 and Y.ndim == 1:
        # Regular grid: use RegularGridInterpolator
        interp_vel = RegularGridInterpolator((Y, X), U, method='linear')
        sample_points = np.column_stack((cN_U, cE_U))
        c_U = interp_vel(sample_points)
    else:
        # Irregular grid: use griddata
        grid_points = np.column_stack((X.ravel(), Y.ravel()))
        vel_values = U.ravel()
        c_U = griddata(grid_points, vel_values, (cE_U, cN_U), method='linear')
    
    return cE_U, cN_U, c_U

# Usage example (matches your Matlab call):
cE_U, cN_U, c_U = improfile(X, Y, U, East, North, 100)

Notes #

• The method='linear' argument matches Matlab's default interpolation for improfile. You can change it to 'nearest' or 'cubic' if needed.
• Make sure the shape of U aligns with your grid: if X is (40,) and Y is (40,), U should be (40,40) where U[i,j] corresponds to (X[j], Y[i]).

内容的提问来源于stack exchange，提问作者Gustavo Gonzalez