如何解决 GDAL 无法识别 NetCDF CF-1.0 坐标变量的问题
GDAL（尤其是 gdalinfo 和 QGIS）在读取符合 CF-1.0 规范的 NetCDF 文件时，可能因坐标变量精度不足而忽略 x/y 坐标轴，导致地理参考丢失；根本原因是使用 float32（f4）存储投影坐标时累积浮点误差破坏了“规则网格”判定逻辑。
gdal（尤其是 `gdalinfo` 和 qgis）在读取符合 cf-1.0 规范的 netcdf 文件时，可能因坐标变量精度不足而忽略 `x`/`y` 坐标轴，导致地理参考丢失；根本原因是使用 `float32`（`f4`）存储投影坐标时累积浮点误差破坏了“规则网格”判定逻辑。
在使用 Python（如 netCDF4 + rasterio）创建带地理参考的 NetCDF 栅格文件时，若严格遵循 CF-1.0 标准定义 projection_x_coordinate（x）和 projection_y_coordinate（y）坐标变量，却在 gdalinfo 中看到角点坐标显示为 (0.0, 0.0) 等非地理值（即未识别为地理坐标），这通常并非元数据缺失，而是 GDAL 对坐标数组数值规律性的隐式校验失败所致。
GDAL 的 NetCDF 驱动在解析 CF 格式时，会检查 x 和 y 变量是否构成严格的一维线性等距网格（即满足 x[i] = x₀ + i × dx，y[j] = y₀ + j × dy）。该判断依赖浮点数的精确表示——当使用 float32（f4）类型存储大范围 UTM 坐标（如 541107.0, 5438536.5）并叠加较小像元尺寸（如 1.5）时，多次累加会产生不可忽略的舍入误差。例如：
# ❌ 危险：float32 累积误差导致 x[1000] ≠ x[0] + 1000 * 1.5（在机器精度内不成立）
x_var = ds.createVariable('x', 'f4', ('x',)) # ← 问题根源
x_var[:] = [541107.0 + i * 1.5 for i in range(7920)] # 实际存储值存在微小偏差
GDAL 检测到 x[1] - x[0] 与 x[2] - x[1] 的差值超出容差（默认约 1e-6 相对误差），便放弃将其视为规则坐标轴，退化为无地理参考的像素坐标系。
✅ 推荐解决方案：升级为 float64
将坐标变量的数据类型从 'f4' 改为 'f8'，可显著提升数值稳定性，确保线性关系在浮点运算中严格保持：
# ✅ 正确：使用 float64 存储坐标，消除累积误差
x_var = ds.createVariable('x', 'f8', ('x',))
y_var = ds.createVariable('y', 'f8', ('y',))
# 生成坐标（仍用 rasterio.transform.xy，但赋值时自动转为 f8）
x_coords = rasterio.transform.xy(transform, cols=list(range(width)), rows=0)[0]
y_coords = rasterio.transform.xy(transform, cols=0, rows=list(range(height)))[1]
x_var[:] = x_coords # 自动以高精度存储
y_var[:] = y_coords
⚠️ 注意事项：
- 即使 xsize/ysize 为整数（如 1.0），若 west/north 原点本身是 float32 表示的十进制数（如 541107.0 在 f4 下实际为 541107.0000000001），仍可能触发误差；
- f8 并非万能：确保 transform 计算本身无精度损失（rasterio.transform.from_origin 返回 Affine 对象，内部使用 float64，安全）；
- CF 属性必须完整：standard_name、units、axis 缺一不可（如 y_var.standard_name = 'projection_y_coordinate'）。
? 备用方案：显式添加 GeoTransform（兼容性更强）
若需兼容旧版 GDAL 或无法修改数据类型，可绕过坐标变量解析，直接在 grid_mapping 变量中写入 GeoTransform 字符串（GDAL 专用扩展属性）：
# 在创建 grid_mapping 变量后添加：
GT0, GT1, GT2, GT3, GT4, GT5 = (
541107.0, # 顶层左像素中心 X
1.5, # 像元宽度（X 方向）
0.0, # 行旋转（通常为 0）
5438536.5, # 顶层左像素中心 Y
0.0, # 列旋转（通常为 0）
-1.5 # 像元高度（Y 方向，负值表示北向上）
)
grid_mapping.GeoTransform = f"{GT0} {GT1} {GT2} {GT3} {GT4} {GT5}"
grid_mapping.spatial_ref = cf_grid_mapping["crs_wkt"] # WKT 定义坐标系
? 提示：GeoTransform 是 GDAL 特有机制，虽非 CF 标准，但被广泛支持；结合 spatial_ref 可确保 QGIS、ArcGIS 等软件正确渲染。
? 总结
GDAL 对 NetCDF CF 坐标识别失败的核心在于浮点精度与规则网格假设的冲突。优先采用 float64 存储坐标变量，是最符合 CF 规范且一劳永逸的修复方式；辅以 GeoTransform 属性则提供强健的后备支持。开发中应始终对地理坐标变量使用双精度，并通过 gdalinfo -stats 验证输出是否显示真实地理角点坐标。