Permalink
Browse files

WCS: support providing creation options that are layer specific

If a Mapfile defines a GEOTIFF outputformat that is enabled in a WCS
layer, then it is possible to define a per-layer value for a creation option,
with a syntax like:

    "wcs_outputformat_GEOTIFF_creationoption_COMPRESS" "LZW"

to say that for 'GEOTIFF' output format, then COMPRESS=LZW should be passed to
the GDAL GTiff driver.

For GDAL 2.3.0 and the GRIB driver that supports band-specific creation options,
it is possible to define such band-specific creation options with for example

   "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)"

In case of band selection, the code will automatically renumber the band number.
  • Loading branch information...
rouault committed Dec 4, 2017
1 parent 834cdb2 commit 6b5ed6e941078a85f743f975d0b90c794807aed0
@@ -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) {
@@ -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 */
@@ -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) {
@@ -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:
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -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
Oops, something went wrong.

0 comments on commit 6b5ed6e

Please sign in to comment.