如何在两个IRIS Cube（GRIB文件）间执行自定义计算？含xarray方案

一、修复纬度坐标不一致问题 #

IRIS对坐标匹配的精度要求极高，哪怕范围一致，坐标点的小数精度、排序方向（如纬度是北→南还是南→北）存在细微差异，都会触发ValueError。可以按以下步骤处理：

• 先排查坐标差异：分别打印两个Cube的纬度坐标，确认数值、精度、顺序是否完全一致：

  print(temp_cube.coord('latitude'))
  print(elev_cube.coord('latitude'))
  

• 强制对齐坐标：直接将海拔Cube的经纬度坐标替换为气温Cube的坐标，确保完全匹配：

  # temp_cube为气温立方体，elev_cube为海拔立方体
  elev_cube.replace_coord(temp_cube.coord('latitude'))
  elev_cube.replace_coord(temp_cube.coord('longitude'))
  

• 若坐标点数量或位置确实存在差异，用重采样对齐：

  import iris
  elev_cube_regridded = elev_cube.regrid(temp_cube, iris.analysis.Linear())
  

二、实现采样点气温计算 #

坐标修复后，按公式批量计算采样点气温：

1. 基础循环实现（适合少量采样点）：

   import numpy as np
   
   # 示例采样点格式：[(经度1, 纬度1, 采样点海拔1), ...]
   sample_points = [(116.3, 39.9, 50), (120.1, 30.2, 100)]
   
   t_sample = np.zeros(len(sample_points))
   for i, (lon, lat, z_sp) in enumerate(sample_points):
       # 插值获取对应点气温（保留时间维度）
       t_air = temp_cube.interpolate([('longitude', lon), ('latitude', lat)], iris.analysis.Linear())
       # 插值获取对应点网格海拔
       z_gb = elev_cube_regridded.interpolate([('longitude', lon), ('latitude', lat)], iris.analysis.Linear())
       # 按公式计算
       t_sample[i] = t_air.data + (z_sp - z_gb.data) * 0.005
   

2. 批量优化实现（适合大量采样点）：

   # 提取采样点的经纬度和海拔数组
   sample_lons = np.array([p[0] for p in sample_points])
   sample_lats = np.array([p[1] for p in sample_points])
   sample_zsp = np.array([p[2] for p in sample_points])
   
   # 批量插值气温（保留72个时间步）
   t_air_all = temp_cube.interpolate([('longitude', sample_lons), ('latitude', sample_lats)], iris.analysis.Linear())
   # 批量插值网格海拔
   z_gb_all = elev_cube_regridded.interpolate([('longitude', sample_lons), ('latitude', sample_lats)], iris.analysis.Linear())
   
   # 广播计算所有时间步的采样点气温
   t_sample_all = t_air_all.data + (sample_zsp[np.newaxis, :] - z_gb_all.data) * 0.005
   

三、是否建议转用xarray？ #

非常建议，尤其适合这类多维度数据的坐标匹配、插值与运算场景：

• xarray对坐标匹配的容错性更高，默认会自动对齐同名坐标，极少出现IRIS式的严格匹配报错；
• 代码更简洁，批量计算的可读性和效率更高，示例代码：

  import xarray as xr
  
  # IRIS Cube转xarray Dataset
  temp_ds = xr.Dataset.from_iris(temp_cube)
  elev_ds = xr.Dataset.from_iris(elev_cube_regridded)
  
  # 采样点转Dataset
  sample_ds = xr.Dataset({
      'z_sp': ('sample', sample_zsp),
      'longitude': ('sample', sample_lons),
      'latitude': ('sample', sample_lats)
  })
  
  # 插值并计算
  t_air_interp = temp_ds['t2m'].interp(longitude=sample_ds['longitude'], latitude=sample_ds['latitude'])
  z_gb_interp = elev_ds['elevation'].interp(longitude=sample_ds['longitude'], latitude=sample_ds['latitude'])
  t_sample = t_air_interp + (sample_ds['z_sp'] - z_gb_interp) * 0.005
  

• xarray生态更完善，配合dask可轻松处理大数据量，可视化、数据导出也更便捷，适合长期的数据分析需求。

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