如何基于地理坐标而非像素坐标提取沿曲线的一维剖面

我拥有地理坐标下的辐射数据（NASA GOLD Mission），每个半球的数据存储在单独文件中。需要针对拼接后的半球数据，沿地理坐标下的10°北磁倾角曲线提取一维剖面。

我尝试参考相关方案，使用scipy.ndimage.map_coordinates实现，但该方法基于像素坐标，导致提取的剖面与地图不匹配；将mapcoords的输出与输入曲线（X）值绘图展示经度变化，黄色剖面也与地图数据不符；使用imshow绘图时，结果也与地图不一致。

相关数据与Python代码已整理好。

原尝试代码：

from netCDF4 import Dataset  
import glob, matplotlib,os;matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy.ma as ma
# 原代码缺失的导入
import numpy as np
import pandas as pd
import scipy.ndimage

Mag = pd.read_csv(r"\mag_cords.csv")
X_dip10s,Y_dip10s = Mag.LON10S.to_numpy(),Mag.LAT10S.to_numpy()
g = Dataset('GOLD_L1C_CH*23*.nc','r'); 
wavelength = g.variables['WAVELENGTH'][:]
radiance = g.variables['RADIANCE'][:]
lat = g.variables['REFERENCE_POINT_LAT'][:] 
lon = g.variables['REFERENCE_POINT_LON'][:] 
g.close()
O5s_ids = np.argwhere((134.3 <= wavelength[50,25,:]) & (wavelength[50,25,:] <= 137.7))
O5s = np.array(np.nansum(radiance[:,:,O5s_ids],axis = 2))*.04 # integrate under the peak!
O5s = np.transpose(O5s[:,:,0])
x=lon[9:-8,9:-8];y=lat[9:-8,9:-8];z=O5s[9:-8,9:-8].T;z=abs(z)#To remove Nans
z2=np.ma.masked_array(z, (z > 30)) 

Zs=scipy.ndimage.map_coordinates(np.transpose(z),
np.vstack((X_dip10s,Y_dip10s)),
mode="nearest")
Xs=scipy.ndimage.map_coordinates(np.transpose(x),
np.vstack((X_dip10s,Y_dip10s)),mode="nearest") 
                                                                                                                         
fig=plt.figure(figsize=(10,10))
axarr = fig.add_subplot(211,projection=ccrs.PlateCarree())
ax2 = fig.add_subplot(212)
grid_lines =axarr.gridlines(draw_labels=True'); axarr.coastlines()
cs1 = axarr.contourf(x,y, z2, transform=ccrs.PlateCarree(),cmap='jet',levels=50)
ax2.plot(X_dip10s,Zs, 'k-', Xs ,Zs, 'k-')

解决方案 #

问题核心是scipy.ndimage.map_coordinates接收的是数组的像素索引坐标，而非地理经纬度坐标。因此必须先把目标地理坐标转换为对应的像素索引，再传入该函数提取剖面。

步骤1：地理坐标转像素索引

假设你的经纬度网格是规则排列的（即x是经度网格，每列对应固定经度；y是纬度网格，每行对应固定纬度），可以用np.interp完成线性插值转换：

# 获取网格的行列数
ny, nx = x.shape

# 提取一维的经度、纬度数组
lon_vals = x[0, :]  # 第一行的经度值，代表所有列的经度
lat_vals = y[:, 0]  # 第一列的纬度值，代表所有行的纬度

# 将目标经纬度转换为像素索引（浮点型，支持插值）
col_indices = np.interp(X_dip10s, lon_vals, np.arange(nx))  # 经度对应列索引
row_indices = np.interp(Y_dip10s, lat_vals, np.arange(ny))  # 纬度对应行索引

# 组合成map_coordinates要求的2xN格式
pixel_coords = np.vstack((row_indices, col_indices))

步骤2：提取剖面数据

用转换后的像素索引调用map_coordinates：

# 提取辐射数据剖面
Zs = scipy.ndimage.map_coordinates(z, pixel_coords, mode="nearest")
# 提取对应位置的经度（用于验证）
Xs = scipy.ndimage.map_coordinates(x, pixel_coords, mode="nearest")

完整修正代码

from netCDF4 import Dataset  
import glob, matplotlib,os;matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy.ma as ma
import numpy as np
import pandas as pd
import scipy.ndimage

Mag = pd.read_csv(r"\mag_cords.csv")
X_dip10s,Y_dip10s = Mag.LON10S.to_numpy(),Mag.LAT10S.to_numpy()
g = Dataset('GOLD_L1C_CH*23*.nc','r'); 
wavelength = g.variables['WAVELENGTH'][:]
radiance = g.variables['RADIANCE'][:]
lat = g.variables['REFERENCE_POINT_LAT'][:] 
lon = g.variables['REFERENCE_POINT_LON'][:] 
g.close()
O5s_ids = np.argwhere((134.3 <= wavelength[50,25,:]) & (wavelength[50,25,:] <= 137.7))
O5s = np.array(np.nansum(radiance[:,:,O5s_ids],axis = 2))*.04 # integrate under the peak!
O5s = np.transpose(O5s[:,:,0])
x=lon[9:-8,9:-8];y=lat[9:-8,9:-8];z=O5s[9:-8,9:-8].T;z=abs(z)#To remove Nans
z2=np.ma.masked_array(z, (z > 30)) 

# ---------------------- 新增：地理坐标转像素索引 ----------------------
ny, nx = x.shape
lon_vals = x[0, :]
lat_vals = y[:, 0]
col_indices = np.interp(X_dip10s, lon_vals, np.arange(nx))
row_indices = np.interp(Y_dip10s, lat_vals, np.arange(ny))
pixel_coords = np.vstack((row_indices, col_indices))
# ---------------------- 提取剖面数据 ----------------------
Zs = scipy.ndimage.map_coordinates(z, pixel_coords, mode="nearest")
Xs = scipy.ndimage.map_coordinates(x, pixel_coords, mode="nearest") 

# 绘图部分保持不变
fig=plt.figure(figsize=(10,10))
axarr = fig.add_subplot(211,projection=ccrs.PlateCarree())
ax2 = fig.add_subplot(212)
grid_lines =axarr.gridlines(draw_labels=True'); axarr.coastlines()
cs1 = axarr.contourf(x,y, z2, transform=ccrs.PlateCarree(),cmap='jet',levels=50)
# 绘制目标曲线和提取的剖面
axarr.plot(X_dip10s, Y_dip10s, 'y-', transform=ccrs.PlateCarree())  # 新增：绘制10°北磁倾角曲线
ax2.plot(X_dip10s,Zs, 'k-', label='目标经度对应剖面')
ax2.plot(Xs,Zs, 'y-', label='提取位置经度对应剖面')
ax2.legend()
plt.show()

关键说明

• 网格匹配：确保lon_vals和lat_vals与你的数据网格维度对应，如果你的经度是按行排列、纬度按列排列，需要调整索引转换逻辑。
• 插值方式：np.interp适用于规则网格；如果是不规则网格，可改用scipy.interpolate.griddata来计算像素索引。
• 维度验证：调用map_coordinates时，输入的数组维度要和像素坐标的行/列对应（比如z的形状是(ny, nx)，则第一维对应行索引，第二维对应列索引）。

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