[GIS] GDAL crashes with segmentation fault when trying to open a NetCDF subdataset with time dimension

gdalnetcdfxarray

I am trying to create a NetCDF file which is compliant with the NetCDF Conventions CF-1.4.

The dimensions of the file are (TIME, Y, X) where Y and X are the polar stereographic grid.

The file that I generate contains three variables: runoff (TIME, Y, X), lat (Y, X) and lon (Y, X). With gdalinfo I can access lat and lon, but not runoff:

$ gdalinfo NETCDF:"myfile.nc":runoff
Segmentation fault (core dumped)

In other words it looks like the TIME dimension is causing problems. I've tried writing all the data in the runoff attribute as zeroes but that didn't help.

The file was created using xarray, with output set to NETCDF4. The ncdump of the file is as follows. Is there anything strange about it that might be causing GDAL to have problems?

netcdf myfile {
dimensions:
    X = 756 ;
    TIME = 12 ;
    Y = 786 ;
variables:
    int64 X(X) ;
        X:units = "meters" ;
        X:standard_name = "projection_x_coordinate" ;
        X:point_spacing = "even" ;
        X:axis = "X" ;
    int64 TIME(TIME) ;
        TIME:standard_name = "time" ;
        TIME:axis = "TIME" ;
        TIME:units = "days since 1958-01-01 00:00:00" ;
        TIME:calendar = "proleptic_gregorian" ;
    int64 Y(Y) ;
        Y:units = "meters" ;
        Y:standard_name = "projection_y_coordinate" ;
        Y:point_spacing = "even" ;
        Y:axis = "Y" ;
    double runoff(TIME, Y, X) ;
        runoff:_FillValue = -9999. ;
        runoff:long_name = "runoff" ;
        runoff:units = "mmWE" ;
        runoff:standard_name = "runoff" ;
        runoff:grid_mapping = "polar_stereographic" ;
    double lon(Y, X) ;
        lon:_FillValue = -9999. ;
        lon:grid_mapping = "polar_stereographic" ;
        lon:units = "degrees_west" ;
        lon:standard_name = "longitude" ;
    byte polar_stereographic ;
        polar_stereographic:grid_mapping_name = "polar_stereographic" ;
        polar_stereographic:scale_factor_at_central_origin = 0. ;
        polar_stereographic:standard_parallel = 70. ;
        polar_stereographic:straight_vertical_longitude_from_pole = -45. ;
        polar_stereographic:latitude_of_projection_origin = 90. ;
        polar_stereographic:false_easting = 0. ;
        polar_stereographic:false_northing = 0. ;
    double lat(Y, X) ;
        lat:_FillValue = -9999. ;
        lat:grid_mapping = "polar_stereographic" ;
        lat:units = "degrees_north" ;
        lat:standard_name = "latitude" ;

// global attributes:
        :Conventions = "CF-1.4" ;
        :history = "redacted" ;
        :institution = "redacted" ;
        :title = "redacted" ;
        :source = "redacted" ;
        :proj4 = "+init=epsg:3413" ;
        :nx = 756. ;
        :ny = 786. ;
        :xmin = -1790000. ;
        :ymax = -65000. ;
        :spacing = 5000.

Best Answer

The problem turned out to be a bug in GDAL, which is now fixed: https://trac.osgeo.org/gdal/ticket/6870.