使用GDAL提取NetCDF多维数组深度层级失败问题排查

背景 #

尝试从NetCDF多维数组文件中提取指定深度层级，但使用gdalmdimtranslate工具时始终失败，报错ERROR 1: Subset specification results in an empty set。查阅GDAL的NetCDF驱动文档后，看到以下关键说明：

维度
NetCDF驱动假设数据遵循UNIDATA的CF-1规范，NetCDF文件内的维度遵循规则：(Z,Y,X)。若维度数量超过3个，驱动会将其合并为波段。例如，对于(P,T,Y,X)类型的4维数组，驱动会将前2个维度相乘（P*T）。驱动会按先递增T再递增P的顺序显示波段，每个波段会附带对应T和P值的元数据。

目标NC文件包含4维数组（变量值、深度、经度、纬度），怀疑问题和上述维度合并机制有关。

执行的命令及错误 #

执行以下两条命令均报错：

gdalmdimtranslate myncfile.nc out.nc -subset 'depth(5727.917)' -array talk
gdalmdimtranslate myncfile.nc out.nc -subset 'depth(5727.91699)' -array talk

其中5727.917和5727.91699是分别通过gdalmdiminfo和ncdump获取的实际深度值。

相关工具输出 #

gdalmdiminfo输出 #

{
  "type": "array",
  "name": "depth",
  "datatype": "Float32",
  "dimensions": [
    {
      "name": "depth",
      "full_name": "/depth",
      "size": 50,
      "type": "VERTICAL",
      "direction": "DOWN",
      "indexing_variable": "/depth"
    }
  ],
  "dimension_size": [
    50
  ],
  "attributes": {
    "units": {
      "datatype": "String",
      "value": "m"
    },
    "positive": {
      "datatype": "String",
      "value": "down"
    },
    "unit_long": {
      "datatype": "String",
      "value": "Meters"
    },
    "long_name": {
      "datatype": "String",
      "value": "Depth"
    },
    "standard_name": {
      "datatype": "String",
      "value": "depth"
    },
    "axis": {
      "datatype": "String",
      "value": "Z"
    },
    "_ChunkSizes": {
      "datatype": "Int32",
      "value": 50
    },
    "_CoordinateAxisType": {
      "datatype": "String",
      "value": "Height"
    },
    "_CoordinateZisPositive": {
      "datatype": "String",
      "value": "down"
    },
    "valid_min": {
      "datatype": "Float32",
      "value": 0.494024992
    },
    "valid_max": {
      "datatype": "Float32",
      "value": 5727.91699
    }
  },
  "unit": "m",
  "values": [0.494024992, 1.54137504, 2.64566898, 3.81949496, 5.07822418, 6.44061422, 7.92956018, 9.57299709, 11.4049997, 13.4671402, 15.81007, 18.4955597, 21.5988197, 25.2114105, 29.4447308, 34.4341507, 40.3440514, 47.3736916, 55.7642899, 65.8072662, 77.8538513, 92.3260727, 109.729301, 130.666, 155.850693, 186.125595, 222.475204, 266.040314, 318.127411, 380.213013, 453.937714, 541.088928, 643.566772, 763.33313, 902.339294, 1062.43994, 1245.29102, 1452.25098, 1684.28406, 1941.89294, 2225.07788, 2533.33594, 2865.70288, 3220.82007, 3597.03198, 3992.48389, 4405.22412, 4833.29102, 5274.78418, 5727.91699]
}

ncdump输出 #

...
data:

 depth = 0.494025, 1.541375, 2.645669, 3.819495, 5.078224, 6.440614, 7.92956,
    9.572997, 11.405, 13.46714, 15.81007, 18.49556, 21.59882, 25.21141,
    29.44473, 34.43415, 40.34405, 47.37369, 55.76429, 65.80727, 77.85385,
    92.32607, 109.7293, 130.666, 155.8507, 186.1256, 222.4752, 266.0403,
    318.1274, 380.213, 453.9377, 541.0889, 643.5668, 763.3331, 902.3393,
    1062.44, 1245.291, 1452.251, 1684.284, 1941.893, 2225.078, 2533.336,
    2865.703, 3220.82, 3597.032, 3992.484, 4405.224, 4833.291, 5274.784,
    5727.917 ;
...

环境信息 #

• GDAL版本：Ubuntu下的3.4.3

排查方向 #

• 浮点数精度匹配问题：GDAL在匹配浮点数维度值时可能存在严格的精度校验，即使值看起来接近，也可能因二进制存储差异导致匹配失败。建议尝试使用索引位置替代具体数值，比如depth(49)（depth数组为0索引，第50个元素的索引是49）。
• 维度合并导致的层级映射问题：由于文件是4维数组，GDAL会将前两个维度合并为波段。需要确认talk变量的维度顺序是否为(depth, time, y, x)或其他组合——如果depth不是被识别为独立的Z维度而是被合并进了波段，那么-subset参数无法直接指定depth值。可以用gdalinfo myncfile.nc查看talk变量的波段元数据，确认每个波段对应的depth值，再通过-b参数指定波段提取。
• CF规范兼容性问题：检查文件是否严格遵循CF-1规范，比如depth变量的axis属性是否正确标记为Z，以及维度顺序是否符合GDAL预期的(Z,Y,X)（或合并后的波段维度）。若维度顺序不符合，GDAL可能无法正确识别depth为垂直维度。
• GDAL版本限制：GDAL 3.4.3可能存在NetCDF子集提取的已知bug，建议尝试升级到较新版本（如3.6+）后重新测试。

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