Merge pull request #5530 from rouault/wcs_per_layer_creation_option

WCS: support providing creation options that are layer specific

mapwcs.c

@@ -1553,6 +1553,81 @@ static int msWCSGetCoverage_ImageCRSSetup(
   return MS_SUCCESS;
 }
 
+/************************************************************************/
+/*                    msWCSApplyLayerCreationOptions()                  */
+/************************************************************************/
+
+void msWCSApplyLayerCreationOptions(layerObj* lp,
+                                    outputFormatObj* format,
+                                    const char* bandlist)
+
+{
+    const char* pszKey;
+    char szKeyBeginning[256];
+    size_t nKeyBeginningLength;
+    int nBands = 0;
+    char** papszBandNumbers = msStringSplit(bandlist, ' ', &nBands);
+
+    snprintf(szKeyBeginning, sizeof(szKeyBeginning),
+             "wcs_outputformat_%s_creationoption_", format->name);
+    nKeyBeginningLength = strlen(szKeyBeginning);
+
+    pszKey = msFirstKeyFromHashTable( &(lp->metadata) );
+    for( ; pszKey != NULL;
+           pszKey = msNextKeyFromHashTable( &(lp->metadata), pszKey) )
+    {
+        if( strncmp(pszKey, szKeyBeginning, nKeyBeginningLength) == 0 )
+        {
+            const char* pszValue = msLookupHashTable( &(lp->metadata), pszKey);
+            const char* pszGDALKey = pszKey + nKeyBeginningLength;
+            if( EQUALN(pszGDALKey, "BAND_", strlen("BAND_")) )
+            {
+                /* Remap BAND specific creation option to the real output
+                 * band number, given the band subset of the request */
+                int nKeyOriBandNumber = atoi(pszGDALKey + strlen("BAND_"));
+                int nTargetBandNumber = -1;
+                int i;
+                for(i = 0; i < nBands; i++ )
+                {
+                    if( nKeyOriBandNumber == atoi(papszBandNumbers[i]) )
+                    {
+                        nTargetBandNumber = i + 1;
+                        break;
+                    }
+                }
+                if( nTargetBandNumber > 0 )
+                {
+                    char szModKey[256];
+                    const char* pszAfterBand =
+                        strchr(pszGDALKey + strlen("BAND_"), '_');
+                    if( pszAfterBand != NULL )
+                    {
+                        snprintf(szModKey, sizeof(szModKey),
+                                 "BAND_%d%s",
+                                 nTargetBandNumber,
+                                 pszAfterBand);
+                        if( lp->debug >= MS_DEBUGLEVEL_VVV ) {
+                            msDebug("Setting GDAL %s=%s creation option\n",
+                                    szModKey, pszValue);
+                        }
+                        msSetOutputFormatOption(format, szModKey, pszValue);
+                    }
+                }
+            }
+            else
+            {
+                if( lp->debug >= MS_DEBUGLEVEL_VVV ) {
+                    msDebug("Setting GDAL %s=%s creation option\n",
+                            pszGDALKey, pszValue);
+                }
+                msSetOutputFormatOption(format, pszGDALKey, pszValue);
+            }
+        }
+    }
+
+    msFreeCharArray( papszBandNumbers, nBands );
+}
+
 /************************************************************************/
 /*                          msWCSGetCoverage()                          */
 /************************************************************************/
@@ -1952,6 +2027,9 @@ this request. Check wcs/ows_enable_request settings.", "msWCSGetCoverage()", par
   msLayerSetProcessingKey(lp, "BANDS", bandlist);
   snprintf(numbands, sizeof(numbands), "%d", msCountChars(bandlist, ',')+1);
   msSetOutputFormatOption(map->outputformat, "BAND_COUNT", numbands);
+
+  msWCSApplyLayerCreationOptions(lp, map->outputformat, bandlist);
+
   free( bandlist );
 
   if(lp->mask) {

mapwcs.h

@@ -102,6 +102,9 @@ void msWCSSetDefaultBandsRangeSetInfo( wcsParamsObj *params,
                                        coverageMetadataObj *cm,
                                        layerObj *lp );
 const char *msWCSGetRequestParameter(cgiRequestObj *request, char *name);
+void msWCSApplyLayerCreationOptions(layerObj* lp,
+                                    outputFormatObj* format,
+                                    const char* bandlist);
 
 /* -------------------------------------------------------------------- */
 /*      Some WCS 1.1 specific functions from mapwcs11.c                 */

mapwcs20.c

@@ -4486,6 +4486,8 @@ this request. Check wcs/ows_enable_request settings.", "msWCSGetCoverage20()", p
   snprintf(numbands, sizeof(numbands), "%d", msCountChars(bandlist, ',')+1);
   msSetOutputFormatOption(map->outputformat, "BAND_COUNT", numbands);
 
+  msWCSApplyLayerCreationOptions(layer, map->outputformat, bandlist);
+
   /* check for the interpolation */
   /* Defaults to NEAREST */
   if(params->interpolation != NULL) {

msautotest/pymod/mstestlib.py

@@ -67,15 +67,53 @@ def get_mapfile_list( argv ):
 
     return map_files
 
+###############################################################################
+# compare_version()
+# Returns 1 if a > b, 0 if a == b, -1 if a < b
+
+def compare_version(version_a, version_b):
+
+    a = version_a.split('.')
+    b = version_b.split('.')
+    while len(a) < 3:
+        a += [ '0' ]
+    while len(b) < 3:
+        b += [ '0' ]
+    a_x, a_y, a_z = int(a[0]), int(a[1]), int(a[2])
+    b_x, b_y, b_z = int(b[0]), int(b[1]), int(b[2])
+    if a_x > b_x:
+        return 1
+    if a_x < b_x:
+        return -1
+    if a_y > b_y:
+        return 1
+    if a_y < b_y:
+        return -1
+    if a_z > b_z:
+        return 1
+    if a_z < b_z:
+        return -1
+    return 0
+
 ###############################################################################
 # has_requires()
 
-def has_requires( version_info, requires_list ):
+def has_requires( version_info, gdal_version, requires_list ):
 
     for item in requires_list:
-        if version_info.find( item ) == -1:
+        if item.startswith('GDAL>='):
+            if gdal_version is None:
+                return 0
+            if compare_version(gdal_version, item[len('GDAL>='):]) < 0:
+                return 0
+        elif item.startswith('GDAL=='):
+            if gdal_version is None:
+                return 0
+            if compare_version(gdal_version, item[len('GDAL=='):]) != 0:
+                return 0
+        elif version_info.find( item ) == -1:
             return 0
-        
+
     return 1
 
 ###############################################################################
@@ -408,6 +446,28 @@ def crlf( filename ):
         else:
             f.write(newdata)
         f.close()
+
+###############################################################################
+def get_gdal_version():
+
+    # First try with GDAL Python bindings, otherwise with gdalinfo binary
+    try:
+        from osgeo import gdal
+        gdal_version = gdal.VersionInfo('VERSION_INFO')
+    except:
+        gdal_version = os.popen( 'gdalinfo --version').read()
+
+    # Parse something like "GDAL x.y.zdev, released..." to extract "x.y.z"
+    if gdal_version.startswith('GDAL '):
+        gdal_version = gdal_version[len('GDAL '):]
+        pos = gdal_version.find('dev')
+        if pos >= 0:
+            return gdal_version[0:pos]
+        pos = gdal_version.find(',')
+        if pos >= 0:
+            return gdal_version[0:pos]
+    return None
+
 ###############################################################################
 # run_tests()
 
@@ -473,12 +533,15 @@ def run_tests( argv ):
     # Get version info.
     version_info = os.popen( shp2img + ' -v' ).read()
     print('version = %s' % version_info)
-
+
+    gdal_version = get_gdal_version()
+    #print('GDAL version = %s' % gdal_version)
+
     ###########################################################################
     # Check directory wide requirements.
     try:
         (runparms_list, requires_list) = read_test_directives( 'all_require.txt' )
-        if not has_requires( version_info, requires_list ):
+        if not has_requires( version_info, gdal_version, requires_list ):
             print('Some or all of the following requirements for this directory of tests\nare not available:')
             print(requires_list)
             return
@@ -511,7 +574,7 @@ def run_tests( argv ):
                 else:
                    runparms_list[i] = ("%s.%s%s"%(resultbase,renderer,resultext),runparms_list[i][1])
 
-        if not has_requires( version_info, requires_list ):
+        if not has_requires( version_info, gdal_version, requires_list ):
             if not quiet:
                 print('    missing some or all of required components, skip.')
             else:

msautotest/wxs/wcs_grib_output.map

@@ -0,0 +1,131 @@
+#
+# Test WCS with GRIB2 output
+#
+# REQUIRES: INPUT=GDAL OUTPUT=PNG SUPPORTS=WCS GDAL>=2.3.0
+
+# RUN_PARMS: wcs_grib_output_wcs10_get_coverage.grb2 [MAPSERV] QUERY_STRING="map=[MAPFILE]&SERVICE=WCS&REQUEST=GetCoverage&VERSION=1.0.0&COVERAGE=temperatures&CRS=EPSG:4326&BBOX=0,-90,180,90&RESX=2.4&RESY=2.4&FORMAT=GRIB" > [RESULT_DEMIME]
+
+# RUN_PARMS: wcs_grib_output_wcs11_get_coverage.grb2 [MAPSERV] QUERY_STRING="map=[MAPFILE]&SERVICE=WCS&REQUEST=GetCoverage&VERSION=1.1.0&IDENTIFIER=temperatures&BOUNDINGBOX=-90,0,90,180,urn:ogc:def:crs:EPSG::4326&GridCS=urn:ogc:def:crs:OGC::imageCRS&GridType=urn:ogc:def:method:WCS:1.1:2dGridIn2dCrs&GridOrigin=90,0&GridOffsets=2.4,-2.4&FORMAT=GRIB" > [RESULT_DEMIME]
+
+# RUN_PARMS: wcs_grib_output_wcs20_get_coverage.grb2 [MAPSERV] QUERY_STRING="map=[MAPFILE]&SERVICE=WCS&REQUEST=GetCoverage&VERSION=2.0.1&COVERAGEID=temperatures&FORMAT=application/x-grib2&SUBSET=long(0,180)&SUBSETTINGCRS=http://www.opengis.net/def/crs/EPSG/0/4326" > [RESULT_DEMIME]
+
+MAP
+
+NAME TEST
+SIZE 400 300
+EXTENT -180 -90 180 90
+MAXSIZE 5000
+
+IMAGETYPE PNG
+TRANSPARENT OFF
+SHAPEPATH "data"
+
+OUTPUTFORMAT
+  NAME GEOTIFF_FLOAT32
+  DRIVER "GDAL/GTiff"
+  MIMETYPE "image/tiff"
+  IMAGEMODE Float32
+  EXTENSION "tif"
+END
+
+OUTPUTFORMAT
+  NAME GRIB
+  DRIVER "GDAL/GRIB"
+  MIMETYPE "application/x-grib2"
+  IMAGEMODE Float32
+  EXTENSION "grib2"
+  FORMATOPTION "DATA_ENCODING=SIMPLE_PACKING"
+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" "*"
+   "ows_srs"	"EPSG:4326"
+
+   # 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"
+   "wcs_formats" "GEOTIFF_FLOAT32 GRIB"
+  END
+END
+
+LAYER
+  NAME temperatures
+  TYPE raster
+  STATUS ON
+  DUMP TRUE
+  DATA "temperatures.tif"
+
+  PROJECTION
+    "init=epsg:4326"
+  END
+  METADATA
+   "ows_extent" "-180 -90 180 90"
+   "wcs_label" "Test label"
+   "ows_srs"	"EPSG:4326"
+   "wcs_resolution" "2.4 2.4"
+   "wcs_bandcount" "1"
+   "wcs_imagemode" "Float32"
+   "wcs_formats" "GEOTIFF_FLOAT32 GRIB"
+   "wcs_description" "Test description"
+   "wcs_metadatalink_href" "http://www.gdal.org/metadata_test_link.html"
+   "wcs_keywordlist" "test,mapserver"
+   "wcs_abstract" "abstract"
+
+   # GDAL creation options for the GRIB output format
+   "wcs_outputformat_GRIB_creationoption_DISCIPLINE" "0(Meteorological)"
+   "wcs_outputformat_GRIB_creationoption_IDS" "Normally this will never be used given the BAND_1_IDS override"
+   "wcs_outputformat_GRIB_creationoption_BAND_1_IDS" "CENTER=54(Montreal) SUBCENTER=0 MASTER_TABLE=4 LOCAL_TABLE=0 SIGNF_REF_TIME=1(Start_of_Forecast) REF_TIME=2017-01-05T00:00:00Z PROD_STATUS=0(Operational) TYPE=2(Analysis_and_forecast)"
+   "wcs_outputformat_GRIB_creationoption_PDS_PDTN" "0"
+   "wcs_outputformat_GRIB_creationoption_PDS_TEMPLATE_ASSEMBLED_VALUES" "0 0 2 47 47 0 0 1 0 103 0 2 255 -127 -2147483647"
+
+   # WCS 2.0 stuff
+   "wcs_native_format" "application/x-grib2"
+   "wcs_band_names" "Temperature_2m"
+   "Temperature_2m_band_interpretation" "interp"
+   "Temperature_2m_band_uom" "Celcius"
+   "Temperature_2m_band_definition"       "Temperature at 2m above ground"
+   "Temperature_2m_band_description"      "Temperature at 2m above ground"
+   "Temperature_2m_interval" "-100 100"
+   "Temperature_2m_significant_figures" "1"
+
+   # WCS 1.x stuff
+   "wcs_nativeformat" "GRIB"
+   "wcs_rangeset_axes" "bands"
+   "wcs_rangeset_name" "Temperatures"
+   "wcs_rangeset_label" "Bands"
+   "wcs_rangeset_description" "Temperatures"
+
+   # "wcs_rangeset_nullvalue" "-99"
+  END
+END
+
+END