WCS and netCDF input/output: automatically create a 3D netCDF file fo…

…r GDAL/netCDF output metadata if the input dataset is itself a 3D netCDF file (if wcs_outputformat_netCDF_mdi_ not defined)

mapwcs.cpp

@@ -2007,13 +2007,94 @@ void msWCSApplySourceDatasetMetadata(layerObj* lp,
         }
         if( !bWCSMetadataFound )
         {
+            int nBands = 0;
+            char** papszBandNumbers = msStringSplit(bandlist, ',', &nBands);
+
+            std::string osExtraDimName;
+            // Special processing if the input dataset is a 3D one
+            {
+                // Check if extra dimensions are declared on the source dataset,
+                // and if so, if there's just a single one.
+                const char* pszDimExtraWithCurl = GDALGetMetadataItem(hDS, "NETCDF_DIM_EXTRA", nullptr);
+                if( pszDimExtraWithCurl &&
+                    strchr(pszDimExtraWithCurl, ',') == nullptr &&
+                    pszDimExtraWithCurl[0] == '{' &&
+                    pszDimExtraWithCurl[strlen(pszDimExtraWithCurl)-1] == '}')
+                {
+                    osExtraDimName.append(
+                        pszDimExtraWithCurl + 1, strlen(pszDimExtraWithCurl) - 2);
+
+                    // Declare the extra dimension name
+                    msSetOutputFormatOption(format,
+                                            "mdi_default_NETCDF_DIM_EXTRA",
+                                            pszDimExtraWithCurl);
+
+                    // Declare the extra dimension definition: size + data type
+                    const char* pszDimExtraDef = GDALGetMetadataItem(hDS,
+                         ("NETCDF_DIM_" + osExtraDimName + "_DEF").c_str(), nullptr);
+                    if( pszDimExtraDef && pszDimExtraDef[0] == '{' &&
+                        pszDimExtraDef[strlen(pszDimExtraDef)-1] == '}' )
+                    {
+                        const auto tokens = msStringSplit(
+                            std::string(pszDimExtraDef + 1, strlen(pszDimExtraDef) -2).c_str(), ',');
+                        if( tokens.size() == 2 )
+                        {
+                            const auto varType = tokens[1];
+                            msSetOutputFormatOption(format,
+                                ("mdi_default_NETCDF_DIM_" + osExtraDimName + "_DEF").c_str(),
+                                (std::string("{") + CPLSPrintf("%d", nBands) + ',' + varType + '}').c_str());
+                        }
+                    }
+
+                    // Declare the extra dimension values
+                    const char* pszDimExtraValues = GDALGetMetadataItem(hDS,
+                         ("NETCDF_DIM_" + osExtraDimName + "_VALUES").c_str(), nullptr);
+                    if( pszDimExtraValues && pszDimExtraValues[0] == '{' &&
+                        pszDimExtraValues[strlen(pszDimExtraValues)-1] == '}' )
+                    {
+                        const auto tokens = msStringSplit(
+                            std::string(pszDimExtraValues + 1, strlen(pszDimExtraValues) -2).c_str(), ',');
+                        if( static_cast<int>(tokens.size()) == GDALGetRasterCount(hDS) )
+                        {
+                            std::string osValue = "{";
+                            for(int i = 0; i < nBands; i++ )
+                            {
+                                int nSrcBand = atoi(papszBandNumbers[i]);
+                                assert( nSrcBand >= 1 && nSrcBand <= static_cast<int>(tokens.size()) );
+                                if( i > 0 )
+                                    osValue += ',';
+                                osValue += tokens[nSrcBand - 1];
+                            }
+                            osValue += '}';
+
+                            msSetOutputFormatOption(format,
+                                ("mdi_default_NETCDF_DIM_" + osExtraDimName + "_VALUES").c_str(),
+                                osValue.c_str());
+                        }
+                    }
+                    else if ( pszDimExtraValues )
+                    {
+                        // If there's a single value
+                        msSetOutputFormatOption(format,
+                                ("mdi_default_NETCDF_DIM_" + osExtraDimName + "_VALUES").c_str(),
+                                pszDimExtraValues);
+                    }
+
+                }
+            }
+
             {
                 char** papszMD = GDALGetMetadata(hDS, NULL);
                 if( papszMD )
                 {
                     for( char** papszIter = papszMD; *papszIter; ++papszIter )
                     {
-                        if( STARTS_WITH(*papszIter, "NC_GLOBAL#") )
+                        // Copy netCDF global attributes, as well as the ones
+                        // of the extra dimension for 3D netCDF files
+                        if( STARTS_WITH(*papszIter, "NC_GLOBAL#") ||
+                            (!osExtraDimName.empty() &&
+                             STARTS_WITH(*papszIter, osExtraDimName.c_str()) &&
+                             (*papszIter)[osExtraDimName.size()] == '#') )
                         {
                             char* pszKey = nullptr;
                             const char* pszValue = CPLParseNameValue(*papszIter, &pszKey);
@@ -2030,10 +2111,7 @@ void msWCSApplySourceDatasetMetadata(layerObj* lp,
                 }
             }
 
-            int nBands = 0;
-            char** papszBandNumbers = msStringSplit(bandlist, ',', &nBands);
-            int i;
-            for(i = 0; i < nBands; i++ )
+            for(int i = 0; i < nBands; i++ )
             {
                 int nSrcBand = atoi(papszBandNumbers[i]);
                 int nDstBand = i + 1;

msautotest/wxs/wcs_netcdf_3d_input_output.map

@@ -0,0 +1,94 @@
+#
+# Test WCS.
+#
+# REQUIRES: INPUT=GDAL OUTPUT=PNG SUPPORTS=WCS
+
+# RUN_PARMS: wcs_netcdf_3d_input_output.nc [MAPSERV] QUERY_STRING="map=[MAPFILE]&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=multi&FORMAT=application/x-netCDF&RANGESUBSET=2,3" > [RESULT_DEMIME][IGNORE_COMPARISON_RESULT_ON_TRAVIS]
+
+MAP
+
+NAME TEST
+SIZE 105 61
+EXTENT 14.4702712 47.8188382 18.0111282 49.8911432
+
+#CONFIG  "MS_ERRORFILE" "stderr"
+
+OUTPUTFORMAT
+  NAME netCDF
+  DRIVER "GDAL/netCDF"
+  MIMETYPE "application/x-netCDF"
+  IMAGEMODE Float32
+  EXTENSION "nc"
+END
+
+PROJECTION
+  "init=epsg:4326"
+END
+
+WEB
+  METADATA
+   # OWS stuff for server
+   "ows_updatesequence"   "2007-10-30T14:23:38Z"
+   "ows_title"            "First Test Service"
+   "ows_fees"             "NONE"
+   "ows_accessconstraints" "NONE"
+   "ows_abstract"         "Test Abstract"
+   "ows_keywordlist"      "keyword,list"
+   "ows_service_onlineresource" "http://198.202.74.215/cgi-bin/wcs_demo"
+   "ows_contactorganization" "OSGeo"
+   "ows_contactperson"    "Frank Warmerdam"
+   "ows_contactposition"  "Software Developer"
+   "ows_contactvoicetelephone" "(613) 754-2041"
+   "ows_contactfacsimiletelephone" "(613) 754-2041x343"
+   "ows_address" "3594 Foymount Rd"
+   "ows_city" "Eganville"
+   "ows_stateorprovince" "Ontario"
+   "ows_postcode" "K0J 1T0"
+   "ows_country" "Canada"
+   "ows_contactelectronicmailaddress" "warmerdam@pobox.com"
+   "ows_hoursofservice" "0800h - 1600h EST"
+   "ows_contactinstructions" "during hours of service"
+   "ows_role" "staff"
+   "ows_enable_request" "*"
+
+   # OGC:WCS
+   "wcs_label"    "Test Label"
+   "wcs_description" "Test description"
+   "wcs_onlineresource"    "http://devgeo.cciw.ca/cgi-bin/mapserv/ecows"
+   "wcs_metadatalink_href" "http://devgeo.cciw.ca/index.html"
+  END #METADATA
+END #WEB
+
+LAYER
+  NAME multi
+  TYPE raster
+  STATUS ON
+
+  DATA expected/wcs_netcdf_3d_output.nc
+
+  PROJECTION
+    "init=epsg:4326"
+  END
+  METADATA
+   "ows_extent" "14.4702712 47.8188382 18.0111282 49.8911432"
+   "wcs_size" "105 61"
+
+   "wcs_label" "Test label"
+   "wcs_formats" "netCDF"
+   "wcs_native_format" "image/tiff"
+   "wcs_description" "Test description"
+   "wcs_metadatalink_href" "http://www.gdal.org/metadata_test_link.html"
+   "wcs_keywordlist" "test,mapserver"
+   "wcs_abstract" "abstract"
+   "wcs_imagemode" "BYTE"
+
+   "wcs_bandcount" "3"
+   "wcs_rangeset_axes" "bands"
+   "wcs_rangeset_name" "name"
+   "wcs_rangeset_label" "Bands"
+   "wcs_rangeset_description" "description"
+   "wcs_rangeset_nullvalue" "0"
+  END
+END
+
+END