如何在contourf/pcolormesh图上提取旋转数据的指定剖面

我把图像存储在numpy数组中，已经编写函数通过将图像索引坐标(i,j)转换为(x,y)并应用旋转矩阵，实现数据按角度theta旋转，返回旋转后的(X,Y)网格。

我需要在同一坐标系中叠加原图与旋转图，并提取特定的垂直、水平剖面，但目前仅能通过map_coordinates函数以'ij'坐标访问旋转图像，尝试scipy的interp2d方法也无效。希望实现在相同xy坐标下提取两张图的同一剖面。

环境配置与函数定义 #

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import pyplot as plt

def rotate_image(arr, dpi, theta_degrees = 0.0, pivot_point = [0,0]):

  theta_radians = (np.pi/180.0)* theta_degrees
  c = round(np.cos(theta_radians), 3)
  s = round(np.sin(theta_radians), 3)

  rotation_matrix = np.array([[c, -s, 0],
                              [s, c, 0],
                              [0, 0,  1]])

  width, height = arr.shape
  pivot_point_xy = np.array([(25.4 / dpi[0])* pivot_point[0], (25.4/dpi[1])*pivot_point[1]])
  pivot_shift_vector = np.array([[pivot_point_xy[0]],
                                 [pivot_point_xy[1]],
                                 [0]])
  
  x = (25.4 / dpi[0]) * np.array(range(width)) # 像素转换为毫米单位
  y = (25.4 / dpi[1]) * np.array(range(height))# 像素转换为毫米单位
  
  XX , YY = np.meshgrid(x,y)
  ZZ = arr
  coordinates = np.stack([XX,YY,ZZ])
  # 平移到旋转中心点、应用旋转、平移回原坐标系
  coordinates_reshape = np.reshape(coordinates, (3,-1))
  translated_coordinates = coordinates_reshape - pivot_shift_vector
  rotated_coordinates = np.matmul(rotation_matrix, translated_coordinates)
  final_coordinates = rotated_coordinates + pivot_shift_vector
  final_coordinates_reshaped = np.reshape(final_coordinates, (3, width, height))
  
  return final_coordinates_reshaped

示例绘图代码 #

img = np.arange(1,26).reshape((5,5))

rotated_img_0 = rotate_image(img, theta_degrees= 0, dpi =[1,1], pivot_point = [2.5,2.5])
rotated_img_1 = rotate_image(img, theta_degrees= 45, dpi =[1,1], pivot_point = [2.5,2.5])

# 绘图
fig, ax = plt.subplots(2, 1, figsize = (10,20))

ax[0].pcolormesh(*rotated_img_0, vmin=0, vmax=rotated_img_0[2].max())
ax[0].pcolormesh(*rotated_img_1, vmin=0, vmax=rotated_img_1[2].max(), alpha = 0.7)
ax[0].hlines(60, rotated_img_1[0].min(), rotated_img_1[0].max() , color = 'black')

ax[1].contourf(*rotated_img_0, vmin=0, vmax=rotated_img_0[2].max())
ax[1].contourf(*rotated_img_1, vmin=0, vmax=rotated_img_1[2].max(), alpha = 0.7)
ax[1].hlines(60, rotated_img_1[0].min(), rotated_img_1[0].max() , color = 'black')

plt.show()

当前困境 #

尝试适配scipy的interp2d方法但对旋转数据无效，map_coordinates仅能通过'ij'坐标处理非旋转数据，简单的i,j切片也仅适用于非旋转数据。需要实现在相同xy坐标下提取两张图的同一剖面。

核心思路是将旋转后的图像重采样到与原图一致的规整XY网格，实现坐标系统一后即可直接提取剖面。

完整实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# 保留原旋转函数不变
def rotate_image(arr, dpi, theta_degrees = 0.0, pivot_point = [0,0]):
    theta_radians = (np.pi/180.0)* theta_degrees
    c = round(np.cos(theta_radians), 3)
    s = round(np.sin(theta_radians), 3)

    rotation_matrix = np.array([[c, -s, 0],
                                [s, c, 0],
                                [0, 0,  1]])

    width, height = arr.shape
    pivot_point_xy = np.array([(25.4 / dpi[0])* pivot_point[0], (25.4/dpi[1])*pivot_point[1]])
    pivot_shift_vector = np.array([[pivot_point_xy[0]],
                                   [pivot_point_xy[1]],
                                   [0]])
    
    x = (25.4 / dpi[0]) * np.arange(width)
    y = (25.4 / dpi[1]) * np.arange(height)
    
    XX , YY = np.meshgrid(x,y)
    ZZ = arr
    coordinates = np.stack([XX,YY,ZZ])
    
    coordinates_reshape = np.reshape(coordinates, (3,-1))
    translated_coordinates = coordinates_reshape - pivot_shift_vector
    rotated_coordinates = np.matmul(rotation_matrix, translated_coordinates)
    final_coordinates = rotated_coordinates + pivot_shift_vector
    final_coordinates_reshaped = np.reshape(final_coordinates, (3, width, height))
    
    return final_coordinates_reshaped

# 生成测试数据
img = np.arange(1,26).reshape((5,5))
rotated_img_0 = rotate_image(img, theta_degrees=0, dpi=[1,1], pivot_point=[2.5,2.5])
rotated_img_1 = rotate_image(img, theta_degrees=45, dpi=[1,1], pivot_point=[2.5,2.5])

# 1. 创建统一采样网格
# 确定两张图覆盖的XY范围
x_min = min(rotated_img_0[0].min(), rotated_img_1[0].min())
x_max = max(rotated_img_0[0].max(), rotated_img_1[0].max())
y_min = min(rotated_img_0[1].min(), rotated_img_1[1].min())
y_max = max(rotated_img_0[1].max(), rotated_img_1[1].max())

# 生成规整网格（采样密度可调整）
xi, yi = np.meshgrid(np.linspace(x_min, x_max, 100), 
                     np.linspace(y_min, y_max, 100))

# 2. 将旋转图重采样到统一网格
points_rot = np.stack([rotated_img_1[0].flatten(), rotated_img_1[1].flatten()], axis=1)
values_rot = rotated_img_1[2].flatten()
z_rot_resampled = griddata(points_rot, values_rot, (xi, yi), method='linear')

# 原图同步映射到统一网格
points_orig = np.stack([rotated_img_0[0].flatten(), rotated_img_0[1].flatten()], axis=1)
values_orig = rotated_img_0[2].flatten()
z_orig_resampled = griddata(points_orig, values_orig, (xi, yi), method='nearest')

# 3. 提取同一XY剖面（以y=60的水平剖面为例）
target_y = 60
# 找到最接近目标y值的网格行
y_idx = np.argmin(np.abs(yi[:,0] - target_y))
# 提取剖面数据
profile_orig = z_orig_resampled[y_idx, :]
profile_rot = z_rot_resampled[y_idx, :]
x_profile = xi[y_idx, :]

# 绘图验证
fig, (ax1, ax2) = plt.subplots(2,1, figsize=(10,12))

# 叠加原图与旋转图
ax1.pcolormesh(xi, yi, z_orig_resampled, vmin=0, vmax=25)
ax1.pcolormesh(xi, yi, z_rot_resampled, vmin=0, vmax=25, alpha=0.7)
ax1.hlines(target_y, x_min, x_max, color='black', lw=2)
ax1.set_title('叠加图与剖面线')

# 绘制剖面对比
ax2.plot(x_profile, profile_orig, label='原图剖面')
ax2.plot(x_profile, profile_rot, label='旋转图剖面', linestyle='--')
ax2.set_title(f'y={target_y}处的水平剖面')
ax2.legend()
plt.tight_layout()
plt.show()

关键说明 #

• griddata支持非结构化点插值，完美适配旋转后不规则的坐标分布，比interp2d更适合该场景。
• 统一网格的采样密度可通过linspace的点数调整，点数越多精度越高，计算量也会相应增加。
• 剖面提取时，通过匹配目标坐标的网格索引直接读取重采样后的数据，实现同一XY坐标下的剖面对比。

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