问询：基于Lambert投影Xarray批量获取大数据集海冰浓度的高效方法

我下载的海冰浓度（sic）数据采用Lambert方位投影，以xc（东向投影坐标）、yc（北向投影坐标）表示，经纬度为二维数组。现有包含大量纬度、经度、时间的DataFrame，基于该Xarray数据集，请问获取对应海冰浓度值的最快/最高效方法是什么？

我的尝试如下：

初始化数据集（实际为超大型DataFrame） #

df_trajs = pd.DataFrame({
    'latitude': [51.5074, 40.7128, 37.7749],
    'longitude': [0.1278, 74.0060, 122.4194],
    'datetime': ['2022-01-01 12:00:00', '2022-01-02 12:00:00', '2022-01-03 12:00:00']
})

df_trajs['datetime'] = pd.to_datetime(df_trajs['datetime'])

实际代码 #

def function_find_var_on_Lambert_Azimuthal_Grid(lat, lon, datetime):    
    dataset = xarray_of_ice_data     
    dataset = dataset.sel(time = datetime, method='nearest') #选择对应时间
    var='ice_conc' #目标变量  
    idxs = (np.abs(ds.lon - lon) + np.abs(ds.lat - lat)).argmin(dim = ["xc", "yc"])   #寻找最近点索引
    value = dataset[var].isel(idxs).values    
    lat_in_ds = ds.lat.isel(idxs).values
    lon_in_ds = ds.lon.isel(idxs).values
    return value

def append_sic_col(df_trajs):
    df_trajs = df_trajs.copy()
    df_trajs.loc[:, 'sic'] = df_trajs[['latitude', 'longitude', 'datetime']].apply(lambda x: function_find_var_on_Lambert_Azimuthal_Grid(*x), axis=1)
    return df_trajs

df_trajs = append_sic_col(df_trajs)

你的现有代码用apply逐行处理，对于超大型DataFrame来说效率极低——逐行操作是Python级别的循环，完全没利用numpy/xarray的向量化优化能力。以下是两种大幅提升效率的方案：

方案1：Xarray向量化插值（推荐） #

直接把DataFrame转换为xarray数据集，一次性完成所有点的插值查找，彻底避免循环：

import xarray as xr

# 假设xarray_of_ice_data是已加载的海冰浓度数据集
ds = xarray_of_ice_data

# 将轨迹DataFrame转换为xarray Dataset，对齐维度
traj_ds = df_trajs.set_index('datetime').to_xarray()

# 用nearest插值（对应原逻辑的最近邻查找）批量获取所有点的海冰浓度
sic_values = ds['ice_conc'].interp(
    time=traj_ds['datetime'],
    lat=traj_ds['latitude'],
    lon=traj_ds['longitude'],
    method='nearest'
).values

# 把结果赋值回原DataFrame
df_trajs['sic'] = sic_values

核心优势：

• 完全依托xarray/numpy的底层C级优化，速度比逐行apply快10~100倍
• 代码简洁，无需自定义循环函数

方案2：预转投影坐标直接索引 #

既然海冰数据用Lambert方位投影存储，把轨迹的经纬度转换成对应的xc/yc投影坐标后，可直接用一维索引查找，速度比二维经纬度匹配更快：

import pyproj

# 定义投影转换（需替换为你的海冰数据集实际投影参数）
proj = pyproj.Proj(proj='laea', lat_0=90, lon_0=0, ellps='WGS84')

# 批量转换轨迹经纬度为xc/yc
df_trajs['xc'], df_trajs['yc'] = proj(df_trajs['longitude'], df_trajs['latitude'])

# 向量化时间+投影坐标索引
sic_values = ds['ice_conc'].sel(
    time=df_trajs['datetime'],
    xc=df_trajs['xc'],
    yc=df_trajs['yc'],
    method='nearest'
).values

df_trajs['sic'] = sic_values

核心优势：

• 投影坐标是一维索引，查找效率远高于二维经纬度数组的匹配
• 同样是向量化操作，无循环开销

原代码低效的关键原因 #

1. apply(axis=1)本质是逐行遍历，每一行都重复执行时间选择、二维数组计算，重复开销极大
2. 每次计算np.abs(ds.lon - lon) + np.abs(ds.lat - lat)都会生成完整的二维数组，内存和计算成本极高

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