NetCDF files should be reconstructed to preserve original variable dimensions

[Mike Taves]

Using tools like gdalinfo and gdalwarp on NetCDF with more than one extra dimensions (beyond 2D band dimension, such as time) will split each into a band. For instance, see variable t from netcdf-4d.nc, which is in trunk/autotest/gdrivers/data.

First, from ncdump -h:

$ ncdump -h netcdf-4d.nc
netcdf netcdf-4d {
dimensions:
        levelist = 2 ;
        longitude = 10 ;
        latitude = 10 ;
        time = UNLIMITED ; // (4 currently)
variables:
        int levelist(levelist) ;
                levelist:long_name = "pressure_level" ;
                levelist:units = "millibars" ;
        float longitude(longitude) ;
                longitude:standard_name = "longitude" ;
                longitude:long_name = "longitude" ;
                longitude:units = "degrees_east" ;
                longitude:axis = "X" ;
        float latitude(latitude) ;
                latitude:standard_name = "latitude" ;
                latitude:long_name = "latitude" ;
                latitude:units = "degrees_north" ;
                latitude:axis = "Y" ;
        double time(time) ;
                time:standard_name = "time" ;
                time:units = "hours since 1900-01-01 00:00:00" ;
                time:calendar = "standard" ;
        int t(time, levelist, latitude, longitude) ;
                t:_FillValue = -32767 ;

// global attributes:
                :Conventions = "CF-1.5" ;
}

Second, from gdalinfo:

$ gdalinfo NetCDF:netcdf-4d.nc:t
Driver: netCDF/Network Common Data Format
Files: netcdf-4d.nc
Size is 10, 10
Coordinate System is `'
Origin = (-1.250000000000000,91.250000000000000)
Pixel Size = (2.500000000000000,-2.500000000000000)
Metadata:
  latitude#axis=Y
  latitude#long_name=latitude
  latitude#standard_name=latitude
  latitude#units=degrees_north
  levelist#long_name=pressure_level
  levelist#units=millibars
  longitude#axis=X
  longitude#long_name=longitude
  longitude#standard_name=longitude
  longitude#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.5
  NETCDF_DIM_EXTRA={time,levelist}
  NETCDF_DIM_levelist_DEF={2,4}
  NETCDF_DIM_levelist_VALUES={925,1000}
  NETCDF_DIM_time_DEF={4,6}
  NETCDF_DIM_time_VALUES={876576,876582,876588,876594}
  t#_FillValue=-32767
  time#calendar=standard
  time#standard_name=time
  time#units=hours since 1900-01-01 00:00:00
Corner Coordinates:
Upper Left  (      -1.250,      91.250)
Lower Left  (  -1.2500000,  66.2500000)
Upper Right (      23.750,      91.250)
Lower Right (  23.7500000,  66.2500000)
Center      (  11.2500000,  78.7500000)
Band 1 Block=10x1 Type=Int32, ColorInterp=Undefined
  NoData Value=-32767
  Metadata:
    _FillValue=-32767
    NETCDF_DIM_levelist=925
    NETCDF_DIM_time=876576
    NETCDF_VARNAME=t
Band 2 Block=10x1 Type=Int32, ColorInterp=Undefined
  NoData Value=-32767
  Metadata:
    _FillValue=-32767
    NETCDF_DIM_levelist=1000
    NETCDF_DIM_time=876576
    NETCDF_VARNAME=t
Band 3 Block=10x1 Type=Int32, ColorInterp=Undefined
  NoData Value=-32767
  Metadata:
    _FillValue=-32767
    NETCDF_DIM_levelist=925
    NETCDF_DIM_time=876582
    NETCDF_VARNAME=t
Band 4 Block=10x1 Type=Int32, ColorInterp=Undefined
  NoData Value=-32767
  Metadata:
    _FillValue=-32767
    NETCDF_DIM_levelist=1000
    NETCDF_DIM_time=876582
    NETCDF_VARNAME=t
Band 5 Block=10x1 Type=Int32, ColorInterp=Undefined
  NoData Value=-32767
  Metadata:
    _FillValue=-32767
    NETCDF_DIM_levelist=925
    NETCDF_DIM_time=876588
    NETCDF_VARNAME=t
Band 6 Block=10x1 Type=Int32, ColorInterp=Undefined
  NoData Value=-32767
  Metadata:
    _FillValue=-32767
    NETCDF_DIM_levelist=1000
    NETCDF_DIM_time=876588
    NETCDF_VARNAME=t
Band 7 Block=10x1 Type=Int32, ColorInterp=Undefined
  NoData Value=-32767
  Metadata:
    _FillValue=-32767
    NETCDF_DIM_levelist=925
    NETCDF_DIM_time=876594
    NETCDF_VARNAME=t
Band 8 Block=10x1 Type=Int32, ColorInterp=Undefined
  NoData Value=-32767
  Metadata:
    _FillValue=-32767
    NETCDF_DIM_levelist=1000
    NETCDF_DIM_time=876594
    NETCDF_VARNAME=t

You can see that GDAL did the right thing: deconstruct t into a series of 2D rasters or bands. There are 8 combinations of time (4) and levelist (2). The NETCDF_DIM_* tags are useful to understand where the indexes for each band are.

However, when reconstructing the NetCDF file from GDAL, each band is stored back as NetCDF Variables. For example, use gdalwarp to reconstruct a new NetCDF file:

$ gdalwarp -of NetCDF NetCDF:netcdf-4d.nc:t from_gdalwarp.nc
Creating output file that is 10P x 10L.
Processing input file NetCDF:netcdf-4d.nc:t.
Using internal nodata values (e.g. -32767) for image NetCDF:netcdf-d.nc:t.
Copying nodata values from source NetCDF:netcdf-4d.nc:t to destination from_gdalwarp.nc.
0...10...20...30...40...50...60...70...80...90...100 - done.
$ ncdump -h from_gdalwarp.nc
netcdf from_gdalwarp {
dimensions:
        lon = 10 ;
        lat = 10 ;
variables:
        int Band1(lat, lon) ;
                Band1:long_name = "GDAL Band Number 1" ;
                Band1:_FillValue = -32767 ;
        int Band2(lat, lon) ;
                Band2:long_name = "GDAL Band Number 2" ;
                Band2:_FillValue = -32767 ;
        int Band3(lat, lon) ;
                Band3:long_name = "GDAL Band Number 3" ;
                Band3:_FillValue = -32767 ;
        int Band4(lat, lon) ;
                Band4:long_name = "GDAL Band Number 4" ;
                Band4:_FillValue = -32767 ;
        int Band5(lat, lon) ;
                Band5:long_name = "GDAL Band Number 5" ;
                Band5:_FillValue = -32767 ;
        int Band6(lat, lon) ;
                Band6:long_name = "GDAL Band Number 6" ;
                Band6:_FillValue = -32767 ;
        int Band7(lat, lon) ;
                Band7:long_name = "GDAL Band Number 7" ;
                Band7:_FillValue = -32767 ;
        int Band8(lat, lon) ;
                Band8:long_name = "GDAL Band Number 8" ;
                Band8:_FillValue = -32767 ;
        double lat(lat) ;
                lat:standard_name = "latitude" ;
                lat:long_name = "latitude" ;
                lat:units = "degrees_north" ;
        double lon(lon) ;
                lon:standard_name = "longitude" ;
                lon:long_name = "longitude" ;
                lon:units = "degrees_east" ;

// global attributes:
                :Conventions = "CF-1.5" ;
                :GDAL = "GDAL 1.11.1, released 2014/09/24" ;
                :history = "Tue Aug 04 09:23:49 2015: GDAL Create( from_gdalwarp
.nc, ... )" ;
}

There are a few unexpected things with the structure of the reconstructed NetCDF file. The dimensions of the original variable t are converted to similar variables for Band1..Band8. It would be expected to have GDAL create a similar int t(time, levelist, latitude, longitude) ; variable, and put the bands back together as they were found. gdalwarp should only change the size of the latitude and longitude dimensions (if truly reprojecting). Also, GDAL renamed these dimension to lat and lon.

The data for the dimensions were stored in the metadata for each band (shown from gdalinfo), but some extra metadata might be required to (e.g.) understand that the time dimension has type int and has UNLIMITED dimension length. The values for these dimensions can be populated from NETCDF_DIM_* tags for each 2D band.

The above was tested using GDAL 1.11.1, released 2014/09/24.

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.