如何从rlon/rlat坐标系中快速提取对应经纬度的字段值？

核心原理说明 #

CORDEX数据集的rlat/rlon是投影后的平面坐标，每个格点对应的真实WGS84经纬度会预先存储为二维lon/lat场，不需要实时做投影解算，所有匹配方案都可以基于预存坐标或投影参数实现高性能批量处理。

高性能批量匹配方案（比逐点KD-Tree快10~100倍） #

方案1：投影坐标直接转换（最优方案，性能最高） #

只要能拿到数据集的投影参数，就可以跳过坐标场匹配，直接做坐标转换后索引取值，无预处理开销，百万级站点匹配仅需百毫秒级耗时。

• CORDEX所有区域的投影参数都存储在数据的rotated_pole变量属性中，直接读取即可
• 用pyproj将站点WGS84经纬度批量转换为对应投影下的x/y坐标，对应rlon/rlat体系
• rlat/rlon本身是等间距一维数组，直接用np.searchsorted即可拿到最近格点索引

import pyproj
import numpy as np
import xarray as xr

# 读取CORDEX数据
ds = xr.open_dataset("你的CORDEX数据路径.nc")
# 读取投影参数构建转换器
crs = pyproj.CRS(ds.rotated_pole.proj4_params)
proj_transformer = pyproj.Transformer.from_crs("EPSG:4326", crs, always_xy=True)
# 读取rlat/rlon一维数组和温度场
rlon = ds.rlon.values
rlat = ds.rlat.values
temp = ds.tas.values # 维度为 [time, rlat, rlon]

# 批量转换站点坐标，site_lons、site_lats为你的站点经纬度数组
site_x, site_y = proj_transformer.transform(site_lons, site_lats)
# 查找最近格点索引，clip避免超出边界
rlon_idx = np.clip(np.searchsorted(rlon, site_x), 0, len(rlon)-1)
rlat_idx = np.clip(np.searchsorted(rlat, site_y), 0, len(rlat)-1)
# 批量提取温度，输出维度为 [站点数, 时间步长]
site_temp = temp[:, rlat_idx, rlon_idx].T

方案2：预生成查找表（适合固定网格+多次匹配场景） #

如果需要反复匹配不同批次的站点，只需要预处理一次网格生成查找表，后续所有匹配都是O(1)耗时，速度最快。

import numpy as np
import xarray as xr

ds = xr.open_dataset("你的CORDEX数据路径.nc")
lon2d = ds.lon.values
lat2d = ds.lat.values
temp = ds.tas.values

# 根据站点精度设置分桶分辨率，0.01度对应约1km精度
lon_res = 0.01
lat_res = 0.01
lon_min, lon_max = lon2d.min(), lon2d.max()
lat_min, lat_max = lat2d.min(), lat2d.max()

# 生成查找表网格，预处理仅执行一次
lon_bins = np.arange(lon_min, lon_max + lon_res, lon_res)
lat_bins = np.arange(lat_min, lat_max + lat_res, lat_res)
lookup_table = np.full((len(lat_bins), len(lon_bins)), fill_value=-1, dtype=int)
for rlat_idx in range(lon2d.shape[0]):
    for rlon_idx in range(lon2d.shape[1]):
        glon, glat = lon2d[rlat_idx, rlon_idx], lat2d[rlat_idx, rlon_idx]
        bin_lon = np.digitize(glon, lon_bins) - 1
        bin_lat = np.digitize(glat, lat_bins) - 1
        if lookup_table[bin_lat, bin_lon] == -1:
            lookup_table[bin_lat, bin_lon] = rlat_idx * lon2d.shape[1] + rlon_idx

# 批量匹配站点，单次查询十万级站点耗时<10毫秒
valid_mask = (site_lons >= lon_min) & (site_lons <= lon_max) & (site_lats >= lat_min) & (site_lats <= lat_max)
valid_lons, valid_lats = site_lons[valid_mask], site_lats[valid_mask]
bin_lons = np.digitize(valid_lons, lon_bins) - 1
bin_lats = np.digitize(valid_lats, lat_bins) - 1
grid_indices = lookup_table[bin_lats, bin_lons]
rlat_indices = grid_indices // lon2d.shape[1]
rlon_indices = grid_indices % lon2d.shape[1]
site_temp = temp[:, rlat_indices, rlon_indices].T

方案3：向量化cKDTree批量查询（适合单次匹配场景） #

你之前用KD-Tree速度慢大概率是逐点循环查询，改用批量查询+多线程可以把速度提升数十倍，适合只需要做一次匹配的场景。

from scipy.spatial import cKDTree
import numpy as np
import xarray as xr

ds = xr.open_dataset("你的CORDEX数据路径.nc")
lon2d = ds.lon.values
lat2d = ds.lat.values
temp = ds.tas.values

# 展平二维坐标场
grid_coords = np.stack([lon2d.ravel(), lat2d.ravel()], axis=1)
tree = cKDTree(grid_coords)
# 一次性批量查询所有站点，n_jobs=-1调用所有CPU核心
site_coords = np.stack([site_lons, site_lats], axis=1)
distances, indices = tree.query(site_coords, k=1, n_jobs=-1)
rlat_indices = indices // lon2d.shape[1]
rlon_indices = indices % lon2d.shape[1]
site_temp = temp[:, rlat_indices, rlon_indices].T

性能对比参考 #

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