如何将经纬度坐标对设置为GOES卫星xarray数据集主索引

问题说明 #

处理自带x/y扫描角坐标的GOES卫星数据时，常规流程为数据转DataFrame后，再转为xarray DataArray开展运算。现有经纬度计算逻辑仅将经纬度作为非索引附加坐标绑定到数据集，这类随两个维度变化的2D成对坐标无法直接设置为数据集索引，导致拆分出DataArray时经纬度信息不会被保留，系统仍默认以原始x/y作为维度坐标。
原有经纬度计算实现如下：

def calc_latlon(ds):
    
    x = ds.x
    y = ds.y
    goes_imager_projection = ds.goes_imager_projection
    
    x,y = np.meshgrid(x,y)
    
    r_eq = goes_imager_projection.attrs["semi_major_axis"]
    r_pol = goes_imager_projection.attrs["semi_minor_axis"]
    l_0 = goes_imager_projection.attrs["longitude_of_projection_origin"] * (np.pi/180)
    h_sat = goes_imager_projection.attrs["perspective_point_height"]
    H = r_eq + h_sat
    
    a = np.sin(x)**2 + (np.cos(x)**2 * (np.cos(y)**2 + (r_eq**2 / r_pol**2) * np.sin(y)**2))
    b = -2 * H * np.cos(x) * np.cos(y)
    c = H**2 - r_eq**2
    
    r_s = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)
    
    s_x = r_s * np.cos(x) * np.cos(y)
    s_y = -r_s * np.sin(x)
    s_z = r_s * np.cos(x) * np.sin(y)
    
    lat = np.arctan((r_eq**2 / r_pol**2) * (s_z / np.sqrt((H-s_x)**2 +s_y**2))) * (180/np.pi)
    lon = (l_0 - np.arctan(s_y / (H-s_x))) * (180/np.pi)
    
    ds = ds.assign_coords({
        "lat":(["y","x"],lat),
        "lon":(["y","x"],lon)
    })
    ds.lat.attrs["units"] = "degrees_north"
    ds.lon.attrs["units"] = "degrees_east"
    return ds

当前数据集元数据显示x/y仍为参考坐标系维度，经纬度仅为附加坐标；拆解为DataArray后系统仍以x/y作为坐标系，未使用已计算的经纬度坐标。此前尝试修改经纬度输出格式、拆分坐标对设置数据集索引、直接通过坐标对设置索引三种思路，均未实现预期效果。

实现方法 #

xarray的维度索引坐标必须为1D结构，且与对应维度严格一一映射，计算得到的2D经纬度网格本身不符合维度索引的要求，因此直接设置索引的操作无法生效。根据使用场景可选择以下两种可落地的方案：

方案1：保留原始x/y维度，标记经纬度为CF规范地理坐标（推荐，无插值误差） #

不需要修改原始维度结构，仅需在绑定经纬度坐标时补全CF公约要求的元数据，即可让xarray及所有支持CF规范的工具（cartopy、rioxarray、NetCDF客户端等）自动识别经纬度为地理参考坐标，拆分DataArray、经纬度切片等操作都不会丢失坐标信息。修改后的计算函数如下：

import numpy as np

def calc_latlon(ds):
    x = ds.x
    y = ds.y
    goes_proj = ds.goes_imager_projection
    x_mesh, y_mesh = np.meshgrid(x, y)

    # 提取GOES投影参数
    r_eq = goes_proj.attrs["semi_major_axis"]
    r_pol = goes_proj.attrs["semi_minor_axis"]
    lon_origin = goes_proj.attrs["longitude_of_projection_origin"] * (np.pi/180)
    h_sat = goes_proj.attrs["perspective_point_height"]
    H = r_eq + h_sat

    # 静止卫星投影几何计算
    a = np.sin(x_mesh)**2 + (np.cos(x_mesh)**2 * (np.cos(y_mesh)**2 + (r_eq**2 / r_pol**2) * np.sin(y_mesh)**2))
    b = -2 * H * np.cos(x_mesh) * np.cos(y_mesh)
    c = H**2 - r_eq**2
    r_s = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)

    s_x = r_s * np.cos(x_mesh) * np.cos(y_mesh)
    s_y = -r_s * np.sin(x_mesh)
    s_z = r_s * np.cos(x_mesh) * np.sin(y_mesh)

    lat = np.arctan((r_eq**2 / r_pol**2) * (s_z / np.sqrt((H-s_x)**2 + s_y**2))) * (180/np.pi)
    lon = (lon_origin - np.arctan(s_y / (H - s_x))) * (180/np.pi)

    # 绑定坐标并补全CF规范元数据
    ds = ds.assign_coords(
        lat=(["y", "x"], lat, {
            "units": "degrees_north",
            "standard_name": "latitude",
            "axis": "Y"
        }),
        lon=(["y", "x"], lon, {
            "units": "degrees_east",
            "standard_name": "longitude",
            "axis": "X"
        })
    )
    return ds

处理完成后可直接通过ds.sel(lat=slice(..., ...), lon=slice(..., ...))按经纬度范围切片，提取的DataArray会自动携带经纬度坐标信息。

方案2：重采样到规则经纬度网格，将经纬度替换为1D维度索引 #

如果需要将经纬度转为可直接作为索引的1D维度（例如与其他规则经纬度网格数据集做匹配运算），可在计算得到经纬度后，通过插值生成规则网格，再重置维度：

# 先通过方案1的函数得到带经纬度坐标的数据集
ds = calc_latlon(ds)

# 按需求定义目标规则经纬度网格的分辨率和范围
target_lon = np.arange(np.floor(ds.lon.min().values), np.ceil(ds.lon.max().values), 0.02)
target_lat = np.arange(np.floor(ds.lat.min().values), np.ceil(ds.lat.max().values), 0.02)

# 插值到规则经纬度网格
ds_regular = ds.interp(lon=target_lon, lat=target_lat).reset_coords(["x", "y"], drop=True)

重采样后的数据集以lat、lon为1D维度索引，拆分出的DataArray默认使用经纬度作为坐标系。

注意事项 #

• 不要强行将2D经纬度网格设置为维度索引，该操作不符合xarray的底层坐标规则，也是此前尝试的三种方案无法生效的核心原因
• 方案1为GOES等静止卫星数据的通用处理方式，不会引入插值误差，也不会丢失原始投影信息，绝大多数场景优先选择该方案
• 补全CF元数据后存储为NetCDF文件，所有常规GIS、遥感软件可直接识别数据的地理参考，无需额外配置投影参数

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