Missing metadata in netCDF output with WCS

[bokhorst]

Hi, I am using MapServer for a WCS service delivering netCDF data (input data is netCDF as well). However, the metadata is not copied from the input to the output. For example, the rather critical NoDataValue (_FillValue in netCDF) is incorrect in the output: it is set to 0.0f instead of -3.e+33f. (examples are provided below)

When I use gdal_translate for the same kind of query, then the output file contains all necessary metadata. This implies that the bug is not in the GDAL driver but in its implementation in Mapserver.

Digging in the source code reveiled that gdal_translate creates a virtual object "VRTDataset" and explicitly copies the metadata from the input object (a "GDALDataSet") by calling GetMetadata(). MapServer in contrast creates a MEMDataSet object but apparently no metadata is copied to this new object. The GetNoDataValue() function is called from the netCDF CreateCopy() function: debugging shows that no result is returned in the case of a Mapserver (hence the 0.0f) and the correct result is returned for gdal_translate().

Original netCDF:

netcdf tide {
dimensions:
        longitude = 3 ;
        latitude = 5 ;
        time = UNLIMITED ; // (1 currently)
variables:
        float longitude(longitude) ;
                longitude:units = "degrees_east" ;
        float latitude(latitude) ;
                latitude:units = "degrees_north" ;
        float time(time) ;
                time:units = "days since 1970-02-01 00:00:00" ;
                time:time_origin = "02-jan-1970 00:00:00" ;
        float z(time, latitude, longitude) ;
                z:units = "m" ;
                z:_FillValue = -3.e+33f ;
}

Mapserver WCS output:

http://xxx/mapserv?map=tidewcs.map&bbox=-180,-76,-179,-78&resx=0.5&resy=-0.5&format=netCDF&request=getCoverage&version=1.0.0&service=WCS&coverage=tide.z&CRS=EPSG:4326

netcdf test {
dimensions:
        x = 2 ;
        y = 4 ;
variables:
        char GDAL_Geographics ;
                GDAL_Geographics:Northernmost_Northing = -78.25 ;
                GDAL_Geographics:Southernmost_Northing = -75.75 ;
                GDAL_Geographics:Easternmost_Easting = -179.25 ;
                GDAL_Geographics:Westernmost_Easting = -180.25 ;
                GDAL_Geographics:spatial_ref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]" ;
                GDAL_Geographics:GeoTransform = "-180.25 0.5 -0 -78.25 0 0.5 " ;
                GDAL_Geographics:grid_mapping_name = "Geographics Coordinate System" ;
                GDAL_Geographics:long_name = "Grid_latitude" ;
        float Band1(y, x) ;
                Band1:_FillValue = 0.f ;
                Band1:grid_mapping = "GDAL_Geographics" ;
                Band1:long_name = "GDAL Band Number 1" ;

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

Gdal_translate output:

gdal_translate -of netCDF -a_srs EPSG:4326  -projwin -180 -78 -179 -76 netcdf:tide.nc:z test.nc

netcdf test {
dimensions:
        x = 2 ;
        y = 4 ;
variables:
        char GDAL_Geographics ;
                GDAL_Geographics:Northernmost_Northing = -78.25 ;
                GDAL_Geographics:Southernmost_Northing = -76.25 ;
                GDAL_Geographics:Easternmost_Easting = -179.25 ;
                GDAL_Geographics:Westernmost_Easting = -180.25 ;
                GDAL_Geographics:spatial_ref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]" ;
                GDAL_Geographics:GeoTransform = "-180.25 0.5 0 -78.25 0 0.5 " ;
                GDAL_Geographics:grid_mapping_name = "Geographics Coordinate System" ;
                GDAL_Geographics:long_name = "Grid_latitude" ;
        float Band1(y, x) ;
                Band1:_FillValue = -3.e+33f ;
                Band1:grid_mapping = "GDAL_Geographics" ;
                Band1:long_name = "GDAL Band Number 1" ;
                Band1:NETCDF_VARNAME = "z" ;
                Band1:NETCDF_DIMENSION_time = "13514.000000" ;
                Band1:NETCDF_time_units = "days since 1970-02-01 00:00:00" ;

// global attributes:
                :Conventions = "CF-1.0" ;
                :z-units = "m" ;
                :z-_FillValue = "-3.000000e+33" ;
                :time-units = "days since 1970-02-01 00:00:00" ;
                :time-time_origin = "02-jan-1970 00:00:00" ;
                :latitude-units = "degrees_north" ;
                :longitude-units = "degrees_east" ;
}