Search code examples
linuxbashnetcdfgdal

NetCDF multidimensional array issue with gdal

I am trying to extract depth levels from a NetCDF multidimensional array file and it seems that those depths cannot be queried with gdal tools (eg. gdalmdimtranslate). After checking the gdal.org website, I found this section under the NetCDF drivers section (https://gdal.org/drivers/raster/netcdf.html):

Dimension The NetCDF driver assume that data follows the CF-1 convention from UNIDATA The dimensions inside the NetCDF file use the following rules: (Z,Y,X). If there are more than 3 dimensions, the driver will merge them into bands. For example if you have an 4 dimension arrays of the type (P, T, Y, X). The driver will multiply the last 2 dimensions (P*T). The driver will display the bands in the following order. It will first increment T and then P. Metadata will be displayed on each band with its corresponding T and P values.

The *.nc file has several bands with 4 dimension arrays (variable value, depth, longitude, latitude). I am not sure why gdal is not recognizing the depth levels, is it a metadata issue related to the above "multiplication of 2 dimensions" explanation? gdal throws the following error message ERROR 1: Subset specification results in an empty set with:

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

Both depth levels (they should be the same value around 5727) are obtained with gdalmdiminfo and ncdump as follows:

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 ;
...

The sample file for testing is here: https://file.io/n8M7J8BLiVsv

Where could it be the problem? Gdal version 3.4.3 in Ubuntu.

Any pointers are welcomed,

Solution

  • Found it, based on the manual (https://gdal.org/drivers/raster/netcdf.html):

    Example of creation of a Time,Z,Y,X 4D file in Python:

# Create in-memory file with required metadata to define the extra >2D
# dimensions
size_z = 2
size_time = 3
src_ds = gdal.GetDriverByName('MEM').Create('', 4, 3, size_z * size_time)
src_ds.SetMetadataItem('NETCDF_DIM_EXTRA', '{time,Z}')
# 6 is NC_DOUBLE
src_ds.SetMetadataItem('NETCDF_DIM_Z_DEF', f"{{{size_z},6}}")
src_ds.SetMetadataItem('NETCDF_DIM_Z_VALUES', '{1.25,2.50}')
src_ds.SetMetadataItem('Z#axis', 'Z')
src_ds.SetMetadataItem('NETCDF_DIM_time_DEF', f"{{{size_time},6}}")
src_ds.SetMetadataItem('NETCDF_DIM_time_VALUES', '{1,2,3}')
src_ds.SetMetadataItem('time#axis', 'T')
src_ds.SetGeoTransform([2,1,0,49,0,-1])

# Create netCDF file
gdal.GetDriverByName('netCDF').CreateCopy('out.nc', src_ds)

    All decimals are needed.