如何将RTOFS NetCDF投影转换为常规经纬度网格并实现海流插值查询

问题背景 #

需构建插值工具，输入实际经纬度后返回对应最近的海表流数值。使用NOAA RTOFS全球海流预报数据，解压后得到3个nc文件，通过xarray打开后发现数据集采用特殊网格：经度范围74.12_{1019.12（非常规±180），纬度范围-78.64}89.98，无法直接对应实际经纬度查询。

核心解决方案：将特殊网格转换为WGS84标准经纬度 #

RTOFS采用Mercator拉伸投影网格，存储的Longitude/Latitude是投影后的坐标，需转换为WGS84（EPSG:4326）标准经纬度才能实现经纬度查询。以下是具体步骤：

1. 确认投影参数 #

先查看数据集元数据获取投影属性：

print(ds.attrs)
# 重点提取投影类型、中央经线、基准面等参数，RTOFS默认参数为Mercator投影、中央经线0°、WGS84基准面

2. 执行投影转换 #

使用pyproj库完成坐标转换：

import pyproj

# 定义RTOFS原始投影
rtofs_proj = pyproj.Proj(proj='merc', lon_0=0, datum='WGS84', units='km')
# 定义目标WGS84经纬度投影
wgs84_proj = pyproj.Proj(proj='latlong', datum='WGS84')

# 获取原始网格坐标数组
x = ds.Longitude.values
y = ds.Latitude.values

# 转换为标准经纬度
lon, lat = pyproj.transform(rtofs_proj, wgs84_proj, x, y)

# 更新数据集坐标，替换为标准经纬度
ds = ds.assign_coords(lon=('node', lon), lat=('node', lat))
# 可选：移除原始投影坐标变量
ds = ds.drop_vars(['Longitude', 'Latitude'])

注：部分RTOFS版本的网格维度为x/y，需根据实际数据集调整维度名称。

3. 构建经纬度插值查询函数 #

利用xarray的interp方法实现经纬度查询：

def get_sea_current(target_lon, target_lat, ds):
    # 采用最近邻插值获取对应海流数值（u为东向流速，v为北向流速）
    interp_result = ds.interp(lon=target_lon, lat=target_lat, method='nearest')
    return interp_result.u.values, interp_result.v.values

# 示例调用
target_lon = 120.5
target_lat = 31.2
u_flow, v_flow = get_sea_current(target_lon, target_lat, ds)
print(f"东向流速：{u_flow} m/s，北向流速：{v_flow} m/s")

4. 验证转换结果 #

通过可视化确认经纬度分布是否符合预期：

import matplotlib.pyplot as plt

plt.scatter(ds.lon.values, ds.lat.values, s=1)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('RTOFS Grid Converted to WGS84')
plt.show()

补充提示 #

• 若投影参数存疑，可参考NOAA官方RTOFS网格定义文档，或用gdalinfo工具读取nc文件的投影元数据。
• 针对时间序列数据，查询时需额外指定目标时间维度参数。

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