/src/MapServer/src/mapwcs20.cpp

Source
/******************************************************************************
 * $Id$
 *
 * Project:  MapServer
 * Purpose:  OpenGIS Web Coverage Server (WCS) 2.0 implementation.
 * Author:   Stephan Meissl <stephan.meissl@eox.at>
 *           Fabian Schindler <fabian.schindler@eox.at>
 *           and the MapServer team.
 *
 ******************************************************************************
 * Copyright (c) 2010, 2011 EOX IT Services GmbH, Austria
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "mapserver-config.h"
#if defined(USE_WCS_SVR)

#include <assert.h>
#include "mapserver.h"
#include "maperror.h"
#include "mapthread.h"
#include "mapows.h"
#include "mapwcs.h"
#include "mapgdal.h"
#include <cmath>
#include <float.h>
#include "gdal.h"
#include "cpl_port.h"
#include "maptime.h"
#include "mapprimitive.h"
#include "cpl_string.h"
#include <string.h>

#if defined(USE_LIBXML2)

#include <libxml/tree.h>
#include "maplibxml2.h"
#include <libxml/parser.h>

#endif /* defined(USE_LIBXML2) */

#include <string>

/************************************************************************/
/*                          msXMLStripIndentation                       */
/************************************************************************/

static void msXMLStripIndentation(char *ptr) {
  /* Remove spaces between > and < to get properly indented result */
  char *afterLastClosingBracket = NULL;
  if (*ptr == ' ')
    afterLastClosingBracket = ptr;
  while (*ptr != '\0') {
    if (*ptr == '<' && afterLastClosingBracket != NULL) {
      memmove(afterLastClosingBracket, ptr, strlen(ptr) + 1);
      ptr = afterLastClosingBracket;
    } else if (*ptr == '>') {
      afterLastClosingBracket = ptr + 1;
    } else if (*ptr != ' ' && *ptr != '\n')
      afterLastClosingBracket = NULL;
    ptr++;
  }
}

/************************************************************************/
/*                   msStringParseInteger()                             */
/*                                                                      */
/*      Tries to parse a string as a integer value. If not possible     */
/*      the value MS_WCS20_UNBOUNDED is stored as value.                */
/*      If no characters could be parsed, MS_FAILURE is returned. If at */
/*      least some characters could be parsed, MS_DONE is returned and  */
/*      only if all characters could be parsed, MS_SUCCESS is returned. */
/************************************************************************/

static int msStringParseInteger(const char *string, int *dest) {
  char *parse_check;
  *dest = (int)strtol(string, &parse_check, 0);
  if (parse_check == string) {
    return MS_FAILURE;
  } else if (parse_check - strlen(string) != string) {
    return MS_DONE;
  }
  return MS_SUCCESS;
}

/************************************************************************/
/*                   msStringParseDouble()                              */
/*                                                                      */
/*      Tries to parse a string as a double value. If not possible      */
/*      the value 0 is stored as value.                                 */
/*      If no characters could be parsed, MS_FAILURE is returned. If at */
/*      least some characters could be parsed, MS_DONE is returned and  */
/*      only if all characters could be parsed, MS_SUCCESS is returned. */
/************************************************************************/

static int msStringParseDouble(const char *string, double *dest) {
  char *parse_check = NULL;
  *dest = strtod(string, &parse_check);
  if (parse_check == string) {
    return MS_FAILURE;
  } else if (parse_check - strlen(string) != string) {
    return MS_DONE;
  }
  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSParseTimeOrScalar20()                         */
/*                                                                      */
/*      Parses a string, either as a time or a scalar value and         */
/*      writes the output into the timeScalarUnion.                     */
/************************************************************************/

static int msWCSParseTimeOrScalar20(timeScalarUnion *u, const char *string) {
  struct tm time;
  if (string) {
    while (*string == ' ')
      string++;
  }

  if (!string || strlen(string) == 0 || !u) {
    msSetError(MS_WCSERR, "Invalid string", "msWCSParseTimeOrScalar20()");
    return MS_WCS20_ERROR_VALUE;
  }
  /* if the string is equal to "*" it means the value
   *  of the interval is unbounded                    */
  if (EQUAL(string, "*")) {
    u->scalar = MS_WCS20_UNBOUNDED;
    u->unbounded = 1;
    return MS_WCS20_UNDEFINED_VALUE;
  }

  /* if returned a valid value, use it */
  if (msStringParseDouble(string, &(u->scalar)) == MS_SUCCESS) {
    return MS_WCS20_SCALAR_VALUE;
  }
  /* otherwise it might be a time value */
  msTimeInit(&time);
  if (msParseTime(string, &time) == MS_TRUE) {
    u->time = mktime(&time);
    return MS_WCS20_TIME_VALUE;
  }
  /* the value could neither be parsed as a double nor as a time value */
  else {
    msSetError(MS_WCSERR,
               "String %s could not be parsed to a time or scalar value",
               "msWCSParseTimeOrScalar20()", string);
    return MS_WCS20_ERROR_VALUE;
  }
}

/************************************************************************/
/*                   msWCSCreateSubsetObj20()                           */
/*                                                                      */
/*      Creates a new wcs20SubsetObj and initializes it to standard     */
/*      values.                                                         */
/************************************************************************/

static wcs20SubsetObjPtr msWCSCreateSubsetObj20() {
  wcs20SubsetObjPtr subset = (wcs20SubsetObjPtr)malloc(sizeof(wcs20SubsetObj));
  MS_CHECK_ALLOC(subset, sizeof(wcs20SubsetObj), NULL);

  subset->axis = NULL;
  subset->crs = NULL;
  subset->min.scalar = subset->max.scalar = MS_WCS20_UNBOUNDED;
  subset->min.unbounded = subset->max.unbounded = 0;
  subset->operation = MS_WCS20_SLICE;

  return subset;
}

/************************************************************************/
/*                   msWCSFreeSubsetObj20()                             */
/*                                                                      */
/*      Frees a wcs20SubsetObj and releases all linked resources.       */
/************************************************************************/

static void msWCSFreeSubsetObj20(wcs20SubsetObjPtr subset) {
  if (NULL == subset) {
    return;
  }
  msFree(subset->axis);
  msFree(subset->crs);
  msFree(subset);
}

/************************************************************************/
/*                   msWCSCreateAxisObj20()                             */
/*                                                                      */
/*      Creates a new wcs20AxisObj and initializes it to standard       */
/*      values.                                                         */
/************************************************************************/

static wcs20AxisObjPtr msWCSCreateAxisObj20() {
  wcs20AxisObj *axis = (wcs20AxisObjPtr)malloc(sizeof(wcs20AxisObj));
  MS_CHECK_ALLOC(axis, sizeof(wcs20AxisObj), NULL);

  axis->name = NULL;
  axis->size = 0;
  axis->resolution = MS_WCS20_UNBOUNDED;
  axis->scale = MS_WCS20_UNBOUNDED;
  axis->resolutionUOM = NULL;
  axis->subset = NULL;

  return axis;
}

/************************************************************************/
/*                   msWCSFreeAxisObj20()                               */
/*                                                                      */
/*      Frees a wcs20AxisObj and releases all linked resources.         */
/************************************************************************/

static void msWCSFreeAxisObj20(wcs20AxisObjPtr axis) {
  if (NULL == axis) {
    return;
  }

  msFree(axis->name);
  msFree(axis->resolutionUOM);
  msWCSFreeSubsetObj20(axis->subset);
  msFree(axis);
}

/************************************************************************/
/*                   msWCSFindAxis20()                                  */
/*                                                                      */
/*      Helper function to retrieve an axis by the name from a params   */
/*      object.                                                         */
/************************************************************************/

static wcs20AxisObjPtr msWCSFindAxis20(wcs20ParamsObjPtr params,
                                       const char *name) {
  int i = 0;
  for (i = 0; i < params->numaxes; ++i) {
    if (EQUAL(params->axes[i]->name, name)) {
      return params->axes[i];
    }
  }
  return NULL;
}

/************************************************************************/
/*                   msWCSInsertAxisObj20()                             */
/*                                                                      */
/*      Helper function to insert an axis object into the axes list of  */
/*      a params object.                                                */
/************************************************************************/

static void msWCSInsertAxisObj20(wcs20ParamsObjPtr params,
                                 wcs20AxisObjPtr axis) {
  params->numaxes++;
  params->axes = (wcs20AxisObjPtr *)msSmallRealloc(
      params->axes, sizeof(wcs20AxisObjPtr) * (params->numaxes));
  params->axes[params->numaxes - 1] = axis;
}

/************************************************************************/
/*                   msWCSCreateParamsObj20()                           */
/*                                                                      */
/*      Creates a new wcs20ParamsObj and initializes it to standard     */
/*      values.                                                         */
/************************************************************************/

wcs20ParamsObjPtr msWCSCreateParamsObj20() {
  wcs20ParamsObjPtr params = (wcs20ParamsObjPtr)malloc(sizeof(wcs20ParamsObj));
  MS_CHECK_ALLOC(params, sizeof(wcs20ParamsObj), NULL);

  params->version = NULL;
  params->request = NULL;
  params->service = NULL;
  params->accept_versions = NULL;
  params->accept_languages = NULL;
  params->sections = NULL;
  params->updatesequence = NULL;
  params->ids = NULL;
  params->width = 0;
  params->height = 0;
  params->resolutionX = MS_WCS20_UNBOUNDED;
  params->resolutionY = MS_WCS20_UNBOUNDED;
  params->scale = MS_WCS20_UNBOUNDED;
  params->scaleX = MS_WCS20_UNBOUNDED;
  params->scaleY = MS_WCS20_UNBOUNDED;
  params->resolutionUnits = NULL;
  params->numaxes = 0;
  params->axes = NULL;
  params->format = NULL;
  params->multipart = 0;
  params->interpolation = NULL;
  params->outputcrs = NULL;
  params->subsetcrs = NULL;
  params->bbox.minx = params->bbox.miny = -DBL_MAX;
  params->bbox.maxx = params->bbox.maxy = DBL_MAX;
  params->range_subset = NULL;
  params->format_options = NULL;

  return params;
}

/************************************************************************/
/*                   msWCSFreeParamsObj20()                             */
/*                                                                      */
/*      Frees a wcs20ParamsObj and releases all linked resources.       */
/************************************************************************/

void msWCSFreeParamsObj20(wcs20ParamsObjPtr params) {
  if (NULL == params) {
    return;
  }

  msFree(params->version);
  msFree(params->request);
  msFree(params->service);
  CSLDestroy(params->accept_versions);
  CSLDestroy(params->accept_languages);
  CSLDestroy(params->sections);
  msFree(params->updatesequence);
  CSLDestroy(params->ids);
  msFree(params->resolutionUnits);
  msFree(params->format);
  msFree(params->interpolation);
  msFree(params->outputcrs);
  msFree(params->subsetcrs);
  while (params->numaxes > 0) {
    params->numaxes -= 1;
    msWCSFreeAxisObj20(params->axes[params->numaxes]);
  }
  msFree(params->axes);
  CSLDestroy(params->range_subset);
  CSLDestroy(params->format_options);
  msFree(params);
}

/************************************************************************/
/*                   msWCSParseSubset20()                               */
/*                                                                      */
/*      Parses several string parameters and fills them into the        */
/*      subset object.                                                  */
/************************************************************************/

static int msWCSParseSubset20(wcs20SubsetObjPtr subset, const char *axis,
                              const char *crs, const char *min,
                              const char *max) {
  int ts1, ts2;
  ts1 = ts2 = MS_WCS20_UNDEFINED_VALUE;

  if (subset == NULL) {
    return MS_FAILURE;
  }

  if (axis == NULL || strlen(axis) == 0) {
    msSetError(MS_WCSERR, "Subset axis is not given.", "msWCSParseSubset20()");
    return MS_FAILURE;
  }

  subset->axis = msStrdup(axis);
  if (crs != NULL) {
    subset->crs = msStrdup(crs);
  }

  /* Parse first (probably only) part of interval/point;
   * check whether its a time value or a scalar value     */
  ts1 = msWCSParseTimeOrScalar20(&(subset->min), min);
  if (ts1 == MS_WCS20_ERROR_VALUE) {
    return MS_FAILURE;
  }

  /* check if its an interval */
  /* if there is a comma, then it is */
  if (max != NULL && strlen(max) > 0) {
    subset->operation = MS_WCS20_TRIM;

    /* Parse the second value of the interval */
    ts2 = msWCSParseTimeOrScalar20(&(subset->max), max);
    if (ts2 == MS_WCS20_ERROR_VALUE) {
      return MS_FAILURE;
    }

    /* if at least one boundary is defined, use that value */
    if ((ts1 == MS_WCS20_UNDEFINED_VALUE) ^ (ts2 == MS_WCS20_UNDEFINED_VALUE)) {
      if (ts1 == MS_WCS20_UNDEFINED_VALUE) {
        ts1 = ts2;
      }
    }
    /* if time and scalar values do not fit, throw an error */
    else if (ts1 != MS_WCS20_UNDEFINED_VALUE &&
             ts2 != MS_WCS20_UNDEFINED_VALUE && ts1 != ts2) {
      msSetError(MS_WCSERR,
                 "Interval error: minimum is a %s value, maximum is a %s value",
                 "msWCSParseSubset20()", ts1 ? "time" : "scalar",
                 ts2 ? "time" : "scalar");
      return MS_FAILURE;
    }
    /* if both min and max are unbounded -> throw an error */
    if (subset->min.unbounded && subset->max.unbounded) {
      msSetError(MS_WCSERR, "Invalid values: no bounds could be parsed",
                 "msWCSParseSubset20()");
      return MS_FAILURE;
    }
  }
  /* there is no second value, therefore it is a point.
   * consequently set the operation to slice */
  else {
    subset->operation = MS_WCS20_SLICE;
    if (ts1 == MS_WCS20_UNDEFINED_VALUE) {
      msSetError(MS_WCSERR, "Invalid point value given",
                 "msWCSParseSubset20()");
      return MS_FAILURE;
    }
  }

  subset->timeOrScalar = ts1;

  /* check whether the min is smaller than the max */
  if (subset->operation == MS_WCS20_TRIM) {
    if (subset->timeOrScalar == MS_WCS20_SCALAR_VALUE &&
        subset->min.scalar == MS_WCS20_UNBOUNDED) {
      subset->min.scalar = -MS_WCS20_UNBOUNDED;
    }

    if (subset->timeOrScalar == MS_WCS20_TIME_VALUE &&
        subset->min.time > subset->max.time) {
      msSetError(MS_WCSERR,
                 "Minimum value of subset axis %s is larger than maximum value",
                 "msWCSParseSubset20()", subset->axis);
      return MS_FAILURE;
    }
    if (subset->timeOrScalar == MS_WCS20_SCALAR_VALUE &&
        subset->min.scalar > subset->max.scalar) {
      msSetError(MS_WCSERR,
                 "Minimum value (%f) of subset axis '%s' is larger than "
                 "maximum value (%f).",
                 "msWCSParseSubset20()", subset->min.scalar, subset->axis,
                 subset->max.scalar);
      return MS_FAILURE;
    }
  }
  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSParseSubsetKVPString20()                      */
/*                                                                      */
/*      Creates a new wcs20SubsetObj, parses a string and fills         */
/*      the parsed values into the object. Returns NULL on failure.     */
/*      Subset string: axis [ , crs ] ( intervalOrPoint )               */
/************************************************************************/

static int msWCSParseSubsetKVPString20(wcs20SubsetObjPtr subset, char *string) {
  char *axis, *crs, *min, *max;

  axis = string;
  crs = NULL;
  min = NULL;
  max = NULL;

  /* find first '(' */
  min = strchr(string, '(');

  /* if min could not be found, the string is invalid */
  if (min == NULL) {
    msSetError(MS_WCSERR, "Invalid axis subset string: '%s'",
               "msWCSParseSubsetKVPString20()", string);
    return MS_FAILURE;
  }
  /* set min to first letter */
  *min = '\0';
  ++min;

  /* cut the trailing ')' */
  if (min[strlen(min) - 1] == ')') {
    min[strlen(min) - 1] = '\0';
  }
  /* look if also a max is defined */
  max = strchr(min, ',');
  if (max != NULL) {
    *max = '\0';
    ++max;
  }

  /* look if also a crs is defined */
  crs = strchr(axis, ',');
  if (crs != NULL) {
    *crs = '\0';
    ++crs;
  }

  return msWCSParseSubset20(subset, axis, crs, min, max);
}

/************************************************************************/
/*                   msWCSParseSizeString20()                           */
/*                                                                      */
/*      Parses a string containing the axis and the size as an integer. */
/*      Size string: axis ( size )                                      */
/************************************************************************/

static int msWCSParseSizeString20(char *string, char *outAxis,
                                  size_t axisStringLen, int *outSize) {
  char *number = NULL;
  char *check = NULL;

  /* find first '(', the character before the number */
  number = strchr(string, '(');

  if (NULL == number) {
    msSetError(MS_WCSERR, "Invalid size parameter value.",
               "msWCSParseSize20()");
    return MS_FAILURE;
  }

  /* cut trailing ')' */
  check = strchr(string, ')');
  if (NULL == check) {
    msSetError(MS_WCSERR, "Invalid size parameter value.",
               "msWCSParseSize20()");
    return MS_FAILURE;
  }
  *number = '\0';
  ++number;
  *check = '\0';

  strlcpy(outAxis, string, axisStringLen);

  /* parse size value */
  if (msStringParseInteger(number, outSize) != MS_SUCCESS) {
    msSetError(MS_WCSERR, "Parameter value '%s' is not a valid integer.",
               "msWCSParseSize20()", number);
    return MS_FAILURE;
  }

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSParseResolutionString20()                     */
/*                                                                      */
/*      Parses a resolution string and returns the axis, the units of   */
/*      measure and the resolution value.                               */
/*      Subset string: axis ( value )                                   */
/************************************************************************/

static int msWCSParseResolutionString20(char *string, char *outAxis,
                                        size_t axisStringLen,
                                        double *outResolution) {
  char *number = NULL;
  char *check = NULL;

  /* find brackets */
  number = strchr(string, '(');

  if (NULL == number) {
    msSetError(MS_WCSERR, "Invalid resolution parameter value : %s.",
               "msWCSParseSize20()", string);
    return MS_FAILURE;
  }

  /* cut trailing ')' */
  check = strchr(string, ')');
  if (NULL == check) {
    msSetError(MS_WCSERR, "Invalid size parameter value.",
               "msWCSParseSize20()");
    return MS_FAILURE;
  }

  *number = '\0';
  ++number;
  *check = '\0';

  strlcpy(outAxis, string, axisStringLen);

  if (msStringParseDouble(number, outResolution) != MS_SUCCESS) {
    *outResolution = MS_WCS20_UNBOUNDED;
    msSetError(MS_WCSERR, "Invalid resolution parameter value : %s.",
               "msWCSParseSize20()", number);
    return MS_FAILURE;
  }

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSParseScaleString20()                          */
/*                                                                      */
/*      Parses a scale string and returns the axis, the and the value.  */
/*      Subset string: axis ( value )                                   */
/************************************************************************/

static int msWCSParseScaleString20(char *string, char *outAxis,
                                   size_t axisStringLen, double *outScale) {
  char *number = NULL;
  char *check = NULL;

  /* find brackets */
  number = strchr(string, '(');

  if (NULL == number) {
    msSetError(MS_WCSERR, "Invalid resolution parameter value : %s.",
               "msWCSParseScaleString20()", string);
    return MS_FAILURE;
  }

  /* cut trailing ')' */
  check = strchr(string, ')');
  if (NULL == check || check < number) {
    msSetError(MS_WCSERR, "Invalid scale parameter value.",
               "msWCSParseScaleString20()");
    return MS_FAILURE;
  }

  *number = '\0';
  ++number;
  *check = '\0';

  strlcpy(outAxis, string, axisStringLen);

  if (msStringParseDouble(number, outScale) != MS_SUCCESS || *outScale <= 0.0) {
    *outScale = MS_WCS20_UNBOUNDED;
    msSetError(MS_WCSERR, "Invalid scale parameter value : %s.",
               "msWCSParseScaleString20()", number);
    return MS_FAILURE;
  }

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSParseResolutionString20()                     */
/*                                                                      */
/*      Parses a resolution string and returns the axis, the units of   */
/*      measure and the resolution value.                               */
/*      Subset string: axis ( value )                                   */
/************************************************************************/

static int msWCSParseScaleExtentString20(char *string, char *outAxis,
                                         size_t axisStringLen, int *outMin,
                                         int *outMax) {
  char *number = NULL;
  char *check = NULL;
  char *colon = NULL;

  /* find brackets */
  number = strchr(string, '(');

  if (NULL == number) {
    msSetError(MS_WCSERR, "Invalid extent parameter value : %s.",
               "msWCSParseScaleExtentString20()", string);
    return MS_FAILURE;
  }

  /* find colon */
  colon = strchr(string, ':');

  if (NULL == colon || colon < number) {
    msSetError(MS_WCSERR, "Invalid extent parameter value : %s.",
               "msWCSParseScaleExtentString20()", string);
    return MS_FAILURE;
  }

  /* cut trailing ')' */
  check = strchr(string, ')');

  if (NULL == check || check < colon) {
    msSetError(MS_WCSERR, "Invalid extent parameter value.",
               "msWCSParseScaleExtentString20()");
    return MS_FAILURE;
  }

  *number = '\0';
  ++number;
  *colon = '\0';
  ++colon;
  *check = '\0';

  strlcpy(outAxis, string, axisStringLen);

  if (msStringParseInteger(number, outMin) != MS_SUCCESS) {
    *outMin = 0;
    msSetError(MS_WCSERR, "Invalid min parameter value : %s.",
               "msWCSParseScaleExtentString20()", number);
    return MS_FAILURE;
  } else if (msStringParseInteger(colon, outMax) != MS_SUCCESS) {
    *outMax = 0;
    msSetError(MS_WCSERR, "Invalid resolution parameter value : %s.",
               "msWCSParseScaleExtentString20()", colon);
    return MS_FAILURE;
  }

  if (*outMin > *outMax) {
    msSetError(MS_WCSERR,
               "Invalid extent: lower part is higher than upper part.",
               "msWCSParseScaleExtentString20()");
    return MS_FAILURE;
  }

  return MS_SUCCESS;
}

#if defined(USE_LIBXML2)
/*
  Utility function to get the first child of a node with a given node name
  */

xmlNodePtr msLibXml2GetFirstChild(xmlNodePtr parent, const char *name) {
  xmlNodePtr node;
  if (!parent || !name) {
    return NULL;
  }

  XML_FOREACH_CHILD(parent, node) {
    XML_LOOP_IGNORE_COMMENT_OR_TEXT(node);
    if (EQUAL((char *)node->name, name)) {
      return node;
    }
  }
  return NULL;
}

/*
  Utility function to get the first child of a node with a given node name and
  namespace.
  */

xmlNodePtr msLibXml2GetFirstChildNs(xmlNodePtr parent, const char *name,
                                    xmlNsPtr ns) {
  xmlNodePtr node;
  if (!parent || !name || !ns) {
    return NULL;
  }

  XML_FOREACH_CHILD(parent, node) {
    XML_LOOP_IGNORE_COMMENT_OR_TEXT(node);
    if (EQUAL((char *)node->name, name) && ns == node->ns) {
      return node;
    }
  }
  return NULL;
}

/*
  Utility function to get the first child of a node with a given node name
  */

xmlNodePtr msLibXml2GetFirstChildElement(xmlNodePtr parent) {
  xmlNodePtr node;
  if (!parent) {
    return NULL;
  }

  XML_FOREACH_CHILD(parent, node) {
    if (node->type == XML_ELEMENT_NODE) {
      return node;
    }
  }
  return NULL;
}
#endif /* defined(USE_LIBXML2) */

/************************************************************************/
/*                   msWCSParseRequest20_XMLGetCapabilities()           */
/*                                                                      */
/*      Parses a DOM element, representing a GetCapabilities-request    */
/*      to a params object.                                             */
/************************************************************************/
#if defined(USE_LIBXML2)
static int msWCSParseRequest20_XMLGetCapabilities(xmlNodePtr root,
                                                  wcs20ParamsObjPtr params) {
  xmlNodePtr child;
  char *content = NULL;
  XML_FOREACH_CHILD(root, child) {
    XML_LOOP_IGNORE_COMMENT_OR_TEXT(child)
    else if (EQUAL((char *)child->name, "AcceptVersions")) {
      xmlNodePtr versionNode = NULL;
      XML_FOREACH_CHILD(child, versionNode) {
        XML_LOOP_IGNORE_COMMENT_OR_TEXT(versionNode);
        XML_ASSERT_NODE_NAME(versionNode, "Version");

        content = (char *)xmlNodeGetContent(versionNode);
        params->accept_versions =
            CSLAddString(params->accept_versions, content);
        xmlFree(content);
      }
    }
    else if (EQUAL((char *)child->name, "Sections")) {
      xmlNodePtr sectionNode = NULL;
      XML_FOREACH_CHILD(child, sectionNode) {
        XML_LOOP_IGNORE_COMMENT_OR_TEXT(sectionNode)
        XML_ASSERT_NODE_NAME(sectionNode, "Section");

        content = (char *)xmlNodeGetContent(sectionNode);
        params->sections = CSLAddString(params->sections, content);
        xmlFree(content);
      }
    }
    else if (EQUAL((char *)child->name, "UpdateSequence")) {
      params->updatesequence = (char *)xmlNodeGetContent(child);
    }
    else if (EQUAL((char *)child->name, "AcceptFormats")) {
      /* Maybe not necessary, since only format is xml.   */
      /* At least ignore it, to not generate an error.    */
    }
    else if (EQUAL((char *)child->name, "AcceptLanguages")) {
      xmlNodePtr languageNode;
      XML_FOREACH_CHILD(child, languageNode) {
        XML_LOOP_IGNORE_COMMENT_OR_TEXT(languageNode)
        XML_ASSERT_NODE_NAME(languageNode, "Language");

        content = (char *)xmlNodeGetContent(languageNode);
        params->accept_languages =
            CSLAddString(params->accept_languages, content);
        xmlFree(content);
      }
    }
    else {
      XML_UNKNOWN_NODE_ERROR(child);
    }
  }
  return MS_SUCCESS;
}
#endif

/************************************************************************/
/*                   msWCSParseRequest20_XMLDescribeCoverage()          */
/*                                                                      */
/*      Parses a DOM element, representing a DescribeCoverage-request   */
/*      to a params object.                                             */
/************************************************************************/
#if defined(USE_LIBXML2)
static int msWCSParseRequest20_XMLDescribeCoverage(xmlNodePtr root,
                                                   wcs20ParamsObjPtr params) {
  xmlNodePtr child;
  char *id;

  XML_FOREACH_CHILD(root, child) {
    XML_LOOP_IGNORE_COMMENT_OR_TEXT(child)
    XML_ASSERT_NODE_NAME(child, "CoverageID");

    /* Node content is the coverage ID */
    id = (char *)xmlNodeGetContent(child);
    if (id == NULL || strlen(id) == 0) {
      msSetError(MS_WCSERR, "CoverageID could not be parsed.",
                 "msWCSParseRequest20_XMLDescribeCoverage()");
      return MS_FAILURE;
    }
    /* insert coverage ID into the list */
    params->ids = CSLAddString(params->ids, (char *)id);
    xmlFree(id);
  }
  return MS_SUCCESS;
}
#endif

/************************************************************************/
/*                   msWCSParseRequest20_XMLGetCoverage()               */
/*                                                                      */
/*      Parses a DOM element, representing a GetCoverage-request to a   */
/*      params object.                                                  */
/************************************************************************/
#if defined(USE_LIBXML2)
static int msWCSParseRequest20_XMLGetCoverage(mapObj *map, xmlNodePtr root,
                                              wcs20ParamsObjPtr params) {
  xmlNodePtr child;
  char *id;

  XML_FOREACH_CHILD(root, child) {
    XML_LOOP_IGNORE_COMMENT_OR_TEXT(child)
    else if (EQUAL((char *)child->name, "CoverageID")) {
      /* Node content is the coverage ID */
      id = (char *)xmlNodeGetContent(child);
      if (id == NULL || strlen(id) == 0) {
        msSetError(MS_WCSERR, "CoverageID could not be parsed.",
                   "msWCSParseRequest20_XMLGetCoverage()");
        return MS_FAILURE;
      }

      /* insert coverage ID into the list */
      params->ids = CSLAddString(params->ids, (char *)id);
      xmlFree(id);
    }
    else if (EQUAL((char *)child->name, "Format")) {
      msFree(params->format);
      params->format = (char *)xmlNodeGetContent(child);
    }
    else if (EQUAL((char *)child->name, "Mediatype")) {
      char *content = (char *)xmlNodeGetContent(child);
      if (content != NULL && (EQUAL(content, "multipart/mixed") ||
                              EQUAL(content, "multipart/related"))) {
        params->multipart = MS_TRUE;
      } else {
        msSetError(MS_WCSERR, "Invalid value '%s' for parameter 'Mediatype'.",
                   "msWCSParseRequest20()", content);
        xmlFree(content);
        return MS_FAILURE;
      }
      xmlFree(content);
    }
    else if (EQUAL((char *)child->name, "DimensionTrim")) {
      wcs20AxisObjPtr axis = NULL;
      wcs20SubsetObjPtr subset = NULL;
      xmlNodePtr node = NULL;
      char *axisName = NULL, *min = NULL, *max = NULL, *crs = NULL;

      /* get strings for axis, min and max */
      XML_FOREACH_CHILD(child, node) {
        XML_LOOP_IGNORE_COMMENT_OR_TEXT(node)
        else if (EQUAL((char *)node->name, "Dimension")) {
          if (axisName != NULL) {
            msSetError(MS_WCSERR, "Parameter 'Dimension' is already set.",
                       "msWCSParseRequest20_XMLGetCoverage()");
            return MS_FAILURE;
          }
          axisName = (char *)xmlNodeGetContent(node);
          crs = (char *)xmlGetProp(node, BAD_CAST "crs");
        }
        else if (EQUAL((char *)node->name, "trimLow")) {
          min = (char *)xmlNodeGetContent(node);
        }
        else if (EQUAL((char *)node->name, "trimHigh")) {
          max = (char *)xmlNodeGetContent(node);
        }
        else {
          msFree(axisName);
          msFree(min);
          msFree(max);
          msFree(crs);
          XML_UNKNOWN_NODE_ERROR(node);
        }
      }
      if (NULL == (subset = msWCSCreateSubsetObj20())) {
        msFree(axisName);
        msFree(min);
        msFree(max);
        msFree(crs);
        return MS_FAILURE;
      }

      /* min and max have to have a value */
      if (min == NULL) {
        min = msStrdup("*");
      }
      if (max == NULL) {
        max = msStrdup("*");
      }
      if (msWCSParseSubset20(subset, axisName, crs, min, max) == MS_FAILURE) {
        msWCSFreeSubsetObj20(subset);
        msWCSException(map, "InvalidSubsetting", "subset", "2.0.1");
        return MS_DONE;
      }

      if (NULL == (axis = msWCSFindAxis20(params, subset->axis))) {
        if (NULL == (axis = msWCSCreateAxisObj20())) {
          msFree(axisName);
          msFree(min);
          msFree(max);
          msFree(crs);
          return MS_FAILURE;
        }
        axis->name = msStrdup(subset->axis);
        msWCSInsertAxisObj20(params, axis);
      }

      axis->subset = subset;

      /* cleanup */
      msFree(axisName);
      msFree(min);
      msFree(max);
      msFree(crs);
    }
    else if (EQUAL((char *)child->name, "DimensionSlice")) {
      msSetError(MS_WCSERR, "Operation '%s' is not supported by MapServer.",
                 "msWCSParseRequest20_XMLGetCoverage()", (char *)child->name);
      return MS_FAILURE;
    }
    else if (EQUAL((char *)child->name, "Size")) {
      wcs20AxisObjPtr axis;
      char *axisName;
      char *content;

      if (NULL ==
          (axisName = (char *)xmlGetProp(child, BAD_CAST "dimension"))) {
        msSetError(MS_WCSERR,
                   "Attribute 'dimension' is missing in element 'Size'.",
                   "msWCSParseRequest20_XMLGetCoverage()");
        return MS_FAILURE;
      }

      if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
        if (NULL == (axis = msWCSCreateAxisObj20())) {
          xmlFree(axisName);
          return MS_FAILURE;
        }
        axis->name = msStrdup(axisName);
        msWCSInsertAxisObj20(params, axis);
      }
      xmlFree(axisName);

      content = (char *)xmlNodeGetContent(child);
      if (msStringParseInteger(content, &(axis->size)) != MS_SUCCESS) {
        xmlFree(content);
        msSetError(MS_WCSERR,
                   "Value of element 'Size' could not "
                   "be parsed to a valid integer.",
                   "msWCSParseRequest20_XMLGetCoverage()");
        return MS_FAILURE;
      }
      xmlFree(content);
    }
    else if (EQUAL((char *)child->name, "Resolution")) {
      wcs20AxisObjPtr axis;
      char *axisName;
      char *content;

      if (NULL ==
          (axisName = (char *)xmlGetProp(child, BAD_CAST "dimension"))) {
        msSetError(MS_WCSERR,
                   "Attribute 'dimension' is missing "
                   "in element 'Resolution'.",
                   "msWCSParseRequest20_XMLGetCoverage()");
        return MS_FAILURE;
      }

      if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
        if (NULL == (axis = msWCSCreateAxisObj20())) {
          xmlFree(axisName);
          return MS_FAILURE;
        }
        axis->name = msStrdup(axisName);
        msWCSInsertAxisObj20(params, axis);
      }
      xmlFree(axisName);

      axis->resolutionUOM = (char *)xmlGetProp(child, BAD_CAST "uom");

      content = (char *)xmlNodeGetContent(child);
      if (msStringParseDouble(content, &(axis->resolution)) != MS_SUCCESS) {
        msSetError(MS_WCSERR,
                   "Value of element 'Resolution' could not "
                   "be parsed to a valid value.",
                   "msWCSParseRequest20_XMLGetCoverage()");
        xmlFree(content);
        return MS_FAILURE;
      }
      xmlFree(content);
    }
    else if (EQUAL((char *)child->name, "Interpolation")) {
      /* Deprecated, use wcs:Extension/int:Interpolation/int:globalInterpolation
       */
      msFree(params->interpolation);
      params->interpolation = (char *)xmlNodeGetContent(child);
    }
    else if (EQUAL((char *)child->name, "OutputCRS")) {
      params->outputcrs = (char *)xmlNodeGetContent(child);
    }
    else if (EQUAL((char *)child->name, "rangeSubset")) {
      /* Deprecated, use wcs:Extension/rsub:RangeSubset */
      xmlNodePtr bandNode = NULL;
      XML_FOREACH_CHILD(child, bandNode) {
        char *content = NULL;
        XML_LOOP_IGNORE_COMMENT_OR_TEXT(bandNode);
        XML_ASSERT_NODE_NAME(bandNode, "band");

        content = (char *)xmlNodeGetContent(bandNode);
        params->range_subset = CSLAddString(params->range_subset, content);
        xmlFree(content);
      }
    }
    else if (EQUAL((char *)child->name, "Extension")) {
      xmlNodePtr extensionNode = NULL;
      XML_FOREACH_CHILD(child, extensionNode) {
        XML_LOOP_IGNORE_COMMENT_OR_TEXT(extensionNode);

        if (EQUAL((char *)extensionNode->name, "Scaling")) {
          xmlNodePtr scaleMethodNode =
              msLibXml2GetFirstChildElement(extensionNode);

          if (EQUAL((char *)scaleMethodNode->name, "ScaleByFactor")) {
            xmlNodePtr scaleFactorNode =
                msLibXml2GetFirstChildElement(scaleMethodNode);
            char *content;
            if (!scaleFactorNode ||
                !EQUAL((char *)scaleFactorNode->name, "scaleFactor")) {
              msSetError(MS_WCSERR, "Missing 'scaleFactor' node.",
                         "msWCSParseRequest20_XMLGetCoverage()");
              return MS_FAILURE;
            }
            content = (char *)xmlNodeGetContent(scaleFactorNode);
            if (msStringParseDouble(content, &(params->scale)) != MS_SUCCESS ||
                params->scale < 0.0) {
              msSetError(MS_WCSERR, "Invalid scaleFactor '%s'.",
                         "msWCSParseRequest20_XMLGetCoverage()", content);
              xmlFree(content);
              return MS_FAILURE;
            }
            xmlFree(content);
          }

          else if (EQUAL((char *)scaleMethodNode->name, "ScaleAxesByFactor")) {
            xmlNodePtr scaleAxisNode, axisNode, scaleFactorNode;
            char *axisName, *content;
            wcs20AxisObjPtr axis;

            XML_FOREACH_CHILD(scaleMethodNode, scaleAxisNode) {
              XML_LOOP_IGNORE_COMMENT_OR_TEXT(scaleAxisNode);

              if (!EQUAL((char *)scaleAxisNode->name, "ScaleAxis")) {
                msSetError(MS_WCSERR, "Invalid ScaleAxesByFactor.",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              /* axis */
              if (NULL ==
                  (axisNode = msLibXml2GetFirstChild(scaleAxisNode, "axis"))) {
                msSetError(MS_WCSERR, "Missing axis node",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }
              axisName = (char *)xmlNodeGetContent(axisNode);
              if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
                if (NULL == (axis = msWCSCreateAxisObj20())) {
                  xmlFree(axisName);
                  return MS_FAILURE;
                }
                axis->name = msStrdup(axisName);
                msWCSInsertAxisObj20(params, axis);
              }
              xmlFree(axisName);

              if (axis->scale != MS_WCS20_UNBOUNDED) {
                msSetError(MS_WCSERR,
                           "scaleFactor was already set for axis '%s'.",
                           "msWCSParseRequest20_XMLGetCoverage()", axis->name);
                return MS_FAILURE;
              }

              /* scaleFactor */
              if (NULL == (scaleFactorNode = msLibXml2GetFirstChild(
                               scaleAxisNode, "scaleFactor"))) {
                msSetError(MS_WCSERR, "Missing scaleFactor node",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              content = (char *)xmlNodeGetContent(scaleFactorNode);
              if (msStringParseDouble(content, &(axis->scale)) != MS_SUCCESS ||
                  axis->scale < 0.0) {
                msSetError(MS_WCSERR, "Invalid scaleFactor '%s'.",
                           "msWCSParseRequest20_XMLGetCoverage()", content);
                xmlFree(content);
                return MS_FAILURE;
              }
              xmlFree(content);
            }
          }

          else if (EQUAL((char *)scaleMethodNode->name, "ScaleToSize")) {
            xmlNodePtr scaleAxisNode, axisNode, targetSizeNode;
            char *axisName, *content;
            wcs20AxisObjPtr axis;

            XML_FOREACH_CHILD(scaleMethodNode, scaleAxisNode) {
              XML_LOOP_IGNORE_COMMENT_OR_TEXT(scaleAxisNode);

              if (!EQUAL((char *)scaleAxisNode->name, "TargetAxisSize")) {
                msSetError(MS_WCSERR, "Invalid ScaleToSize.",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              /* axis */
              if (NULL ==
                  (axisNode = msLibXml2GetFirstChild(scaleAxisNode, "axis"))) {
                msSetError(MS_WCSERR, "Missing axis node",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }
              axisName = (char *)xmlNodeGetContent(axisNode);
              if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
                if (NULL == (axis = msWCSCreateAxisObj20())) {
                  xmlFree(axisName);
                  return MS_FAILURE;
                }
                axis->name = msStrdup(axisName);
                msWCSInsertAxisObj20(params, axis);
              }
              xmlFree(axisName);

              if (axis->size != 0) {
                msSetError(MS_WCSERR,
                           "targetSize was already set for axis '%s'.",
                           "msWCSParseRequest20_XMLGetCoverage()", axis->name);
                return MS_FAILURE;
              }

              /* targetSize */
              if (NULL == (targetSizeNode = msLibXml2GetFirstChild(
                               scaleAxisNode, "targetSize"))) {
                msSetError(MS_WCSERR, "Missing targetSize node",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              content = (char *)xmlNodeGetContent(targetSizeNode);
              if (msStringParseInteger(content, &(axis->size)) != MS_SUCCESS ||
                  axis->size <= 0) {
                msSetError(MS_WCSERR, "Invalid targetSize '%s'.",
                           "msWCSParseRequest20_XMLGetCoverage()", content);
                xmlFree(content);
                return MS_FAILURE;
              }
              xmlFree(content);
            }
          }

          else if (EQUAL((char *)scaleMethodNode->name, "ScaleToExtent")) {
            xmlNodePtr scaleAxisNode, axisNode, lowNode, highNode;
            char *axisName, *content;
            wcs20AxisObjPtr axis;
            int low, high;

            XML_FOREACH_CHILD(scaleMethodNode, scaleAxisNode) {
              XML_LOOP_IGNORE_COMMENT_OR_TEXT(scaleAxisNode);

              if (!EQUAL((char *)scaleAxisNode->name, "TargetAxisExtent")) {
                msSetError(MS_WCSERR, "Invalid ScaleToExtent.",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              /* axis */
              if (NULL ==
                  (axisNode = msLibXml2GetFirstChild(scaleAxisNode, "axis"))) {
                msSetError(MS_WCSERR, "Missing axis node",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }
              axisName = (char *)xmlNodeGetContent(axisNode);
              if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
                if (NULL == (axis = msWCSCreateAxisObj20())) {
                  xmlFree(axisName);
                  return MS_FAILURE;
                }
                axis->name = msStrdup(axisName);
                msWCSInsertAxisObj20(params, axis);
              }
              xmlFree(axisName);

              if (axis->size != 0) {
                msSetError(MS_WCSERR,
                           "targetSize was already set for axis '%s'.",
                           "msWCSParseRequest20_XMLGetCoverage()", axis->name);
                return MS_FAILURE;
              }

              /* targetSize */
              if (NULL ==
                  (lowNode = msLibXml2GetFirstChild(scaleAxisNode, "low"))) {
                msSetError(MS_WCSERR, "Missing low node",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              if (NULL ==
                  (highNode = msLibXml2GetFirstChild(scaleAxisNode, "high"))) {
                msSetError(MS_WCSERR, "Missing high node",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              content = (char *)xmlNodeGetContent(lowNode);
              if (msStringParseInteger(content, &low) != MS_SUCCESS) {
                msSetError(MS_WCSERR, "Invalid low value '%s'.",
                           "msWCSParseRequest20_XMLGetCoverage()", content);
                xmlFree(content);
                return MS_FAILURE;
              }
              xmlFree(content);

              content = (char *)xmlNodeGetContent(highNode);
              if (msStringParseInteger(content, &high) != MS_SUCCESS) {
                msSetError(MS_WCSERR, "Invalid high value '%s'.",
                           "msWCSParseRequest20_XMLGetCoverage()", content);
                xmlFree(content);
                return MS_FAILURE;
              }
              xmlFree(content);

              if (high <= low) {
                msSetError(MS_WCSERR, "Invalid extent, high is lower than low.",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              axis->size = high - low;
            }
          }
        }

        /* Range Subset */
        else if (EQUAL((char *)extensionNode->name, "RangeSubset")) {
          xmlNodePtr rangeItemNode = NULL;

          XML_FOREACH_CHILD(extensionNode, rangeItemNode) {
            xmlNodePtr rangeItemNodeChild =
                msLibXml2GetFirstChildElement(rangeItemNode);
            XML_LOOP_IGNORE_COMMENT_OR_TEXT(rangeItemNode);

            XML_ASSERT_NODE_NAME(rangeItemNode, "RangeItem");

            if (!rangeItemNodeChild) {
              msSetError(MS_WCSERR, "Missing RangeComponent or RangeInterval.",
                         "msWCSParseRequest20_XMLGetCoverage()");
              return MS_FAILURE;
            } else if (EQUAL((char *)rangeItemNodeChild->name,
                             "RangeComponent")) {
              char *content = (char *)xmlNodeGetContent(rangeItemNodeChild);
              params->range_subset =
                  CSLAddString(params->range_subset, content);
              xmlFree(content);
            } else if (EQUAL((char *)rangeItemNodeChild->name,
                             "RangeInterval")) {
              xmlNodePtr intervalNode = rangeItemNodeChild;
              xmlNodePtr startComponentNode =
                  msLibXml2GetFirstChild(intervalNode, "startComponent");
              xmlNodePtr endComponentNode =
                  msLibXml2GetFirstChild(intervalNode, "endComponent");
              char *start;
              char *stop;

              if (!startComponentNode || !endComponentNode) {
                msSetError(MS_WCSERR, "Wrong RangeInterval.",
                           "msWCSParseRequest20_XMLGetCoverage()");
                return MS_FAILURE;
              }

              start = (char *)xmlNodeGetContent(startComponentNode);
              stop = (char *)xmlNodeGetContent(endComponentNode);

              std::string value(start);
              value += ':';
              value += stop;

              xmlFree(start);
              xmlFree(stop);

              params->range_subset =
                  CSLAddString(params->range_subset, value.c_str());
            }
          }
        }

        else if (EQUAL((char *)extensionNode->name, "subsettingCrs")) {
          msFree(params->subsetcrs);
          params->subsetcrs = (char *)xmlNodeGetContent(extensionNode);
        }

        else if (EQUAL((char *)extensionNode->name, "outputCrs")) {
          msFree(params->outputcrs);
          params->outputcrs = (char *)xmlNodeGetContent(extensionNode);
        }

        else if (EQUAL((char *)extensionNode->name, "Interpolation")) {
          xmlNodePtr globalInterpolation =
              msLibXml2GetFirstChild(extensionNode, "globalInterpolation");
          char *content;
          if (globalInterpolation == NULL) {
            msSetError(MS_WCSERR, "Missing 'globalInterpolation' node.",
                       "msWCSParseRequest20_XMLGetCoverage()");
            return MS_FAILURE;
          }
          content = (char *)xmlNodeGetContent(globalInterpolation);
          msFree(params->interpolation);
          /* TODO: use URIs/URLs once they are specified */
          params->interpolation = msStrdup(content);
          xmlFree(content);
        }

        /* GeoTIFF parameters */
        else if (EQUAL((char *)extensionNode->name, "parameters") &&
                 extensionNode->ns &&
                 EQUAL((char *)extensionNode->ns->href,
                       "http://www.opengis.net/gmlcov/geotiff/1.0")) {

          xmlNodePtr parameter;

          XML_FOREACH_CHILD(extensionNode, parameter) {
            char *content;
            XML_LOOP_IGNORE_COMMENT_OR_TEXT(parameter);

            content = (char *)xmlNodeGetContent(parameter);

            params->format_options = CSLAddNameValue(
                params->format_options, (char *)parameter->name, content);
            xmlFree(content);
          }
        }
      }
    }
    else {
      XML_UNKNOWN_NODE_ERROR(child);
    }
  }
  return MS_SUCCESS;
}
#endif

/************************************************************************/
/*                   msWCSParseRequest20()                              */
/*                                                                      */
/*      Parses a CGI-request to a WCS 20 params object. It is           */
/*      either a POST or a GET request. In case of a POST request       */
/*      the xml content has to be parsed to a DOM structure             */
/*      before the parameters can be extracted.                         */
/************************************************************************/

int msWCSParseRequest20(mapObj *map, cgiRequestObj *request,
                        owsRequestObj *ows_request, wcs20ParamsObjPtr params) {
  int i;
  if (params == NULL || request == NULL || ows_request == NULL) {
    msSetError(MS_WCSERR, "Internal error.", "msWCSParseRequest20()");
    return MS_FAILURE;
  }

  /* Copy arbitrary service, version and request. */
  params->service = msStrdup(ows_request->service);
  if (ows_request->version != NULL) {
    params->version = msStrdup(ows_request->version);
  }
  params->request = msStrdup(ows_request->request);

  /* Parse the POST request */
  if (request->type == MS_POST_REQUEST) {
#if defined(USE_LIBXML2)
    xmlDocPtr doc = static_cast<xmlDocPtr>(ows_request->document);
    xmlNodePtr root = NULL;
    const char *validate;
    int ret = MS_SUCCESS;

    /* parse to DOM-Structure and get root element */
    if (doc == NULL) {
      const xmlError *error = xmlGetLastError();
      msSetError(MS_WCSERR, "XML parsing error: %s", "msWCSParseRequest20()",
                 error->message);
      return MS_FAILURE;
    }

    root = xmlDocGetRootElement(doc);

    validate = msOWSLookupMetadata(&(map->web.metadata), "CO", "validate_xml");
    if (validate != NULL && EQUAL(validate, "TRUE")) {
      char *schema_dir = msStrdup(
          msOWSLookupMetadata(&(map->web.metadata), "CO", "schemas_dir"));
      if (schema_dir != NULL &&
          (params->version == NULL || EQUALN(params->version, "2.0", 3))) {
        schema_dir = msStringConcatenate(schema_dir, "wcs/2.0.0/wcsAll.xsd");
        if (msOWSSchemaValidation(schema_dir, request->postrequest) != 0) {
          msSetError(MS_WCSERR,
                     "Invalid POST request. "
                     "XML is not valid",
                     "msWCSParseRequest20()");
          return MS_FAILURE;
        }
      }
      msFree(schema_dir);
    }

    if (EQUAL(params->request, "GetCapabilities")) {
      ret = msWCSParseRequest20_XMLGetCapabilities(root, params);
    } else if (params->version != NULL && EQUALN(params->version, "2.0", 3)) {
      if (EQUAL(params->request, "DescribeCoverage")) {
        ret = msWCSParseRequest20_XMLDescribeCoverage(root, params);
      } else if (EQUAL(params->request, "GetCoverage")) {
        ret = msWCSParseRequest20_XMLGetCoverage(map, root, params);
      }
    }
    return ret;

#else  /* defined(USE_LIBXML2) */
    /* TODO: maybe with CPLXML? */
    return MS_FAILURE;
#endif /* defined(USE_LIBXML2) */
  }

  /* Parse the KVP GET request */
  for (i = 0; i < request->NumParams; ++i) {
    char *key = NULL, *value = NULL;
    char **tokens;
    int num, j;
    key = request->ParamNames[i];
    value = request->ParamValues[i];

    if (EQUAL(key, "VERSION")) {
      continue;
    } else if (EQUAL(key, "REQUEST")) {
      continue;
    } else if (EQUAL(key, "SERVICE")) {
      continue;
    } else if (EQUAL(key, "ACCEPTVERSIONS")) {
      tokens = msStringSplit(value, ',', &num);
      for (j = 0; j < num; ++j) {
        params->accept_versions =
            CSLAddString(params->accept_versions, tokens[j]);
      }
      msFreeCharArray(tokens, num);
    } else if (EQUAL(key, "SECTIONS")) {
      tokens = msStringSplit(value, ',', &num);
      for (j = 0; j < num; ++j) {
        params->sections = CSLAddString(params->sections, tokens[j]);
      }
      msFreeCharArray(tokens, num);
    } else if (EQUAL(key, "UPDATESEQUENCE")) {
      msFree(params->updatesequence);
      params->updatesequence = msStrdup(value);
    } else if (EQUAL(key, "ACCEPTFORMATS")) {
      /* ignore */
    } else if (EQUAL(key, "ACCEPTLANGUAGES")) {
      if (params->accept_languages != NULL) {
        CSLDestroy(params->accept_languages);
      }
      params->accept_languages = CSLTokenizeString2(value, ",", 0);
    } else if (EQUAL(key, "COVERAGEID")) {
      if (params->ids != NULL) {
        msSetError(MS_WCSERR,
                   "Parameter 'CoverageID' is already set. "
                   "For multiple IDs use a comma separated list.",
                   "msWCSParseRequest20()");
        return MS_FAILURE;
      }
      params->ids = CSLTokenizeString2(value, ",", 0);
    } else if (EQUAL(key, "FORMAT")) {
      msFree(params->format);
      params->format = msStrdup(value);
    } else if (EQUAL(key, "MEDIATYPE")) {
      if (EQUAL(value, "multipart/mixed") ||
          EQUAL(value, "multipart/related")) {
        params->multipart = MS_TRUE;
      } else {
        msSetError(MS_WCSERR, "Invalid value '%s' for parameter 'Mediatype'.",
                   "msWCSParseRequest20()", value);
        return MS_FAILURE;
      }
    } else if (EQUAL(key, "INTERPOLATION")) {
      msFree(params->interpolation);
      params->interpolation = msStrdup(value);
    } else if (EQUAL(key, "OUTPUTCRS")) {
      msFree(params->outputcrs);
      params->outputcrs = msStrdup(value);
    } else if (EQUAL(key, "SUBSETTINGCRS")) {
      msFree(params->subsetcrs);
      params->subsetcrs = msStrdup(value);
    } else if (EQUAL(key, "SCALEFACTOR")) {
      double scale = MS_WCS20_UNBOUNDED;
      if (params->scale != MS_WCS20_UNBOUNDED) {
        msSetError(MS_WCSERR, "Parameter 'SCALEFACTOR' already set.",
                   "msWCSParseRequest20()");
        return MS_FAILURE;
      } else if (msStringParseDouble(value, &scale) != MS_SUCCESS) {
        msSetError(MS_WCSERR, "Could not parse parameter 'SCALEFACTOR'.",
                   "msWCSParseRequest20()");
        return MS_FAILURE;
      } else if (scale <= 0.0) {
        msSetError(MS_WCSERR, "Invalid value for 'SCALEFACTOR'.",
                   "msWCSParseRequest20()");
        return MS_FAILURE;
      }
      params->scale = scale;
    } else if (EQUAL(key, "SCALEAXES")) {
      wcs20AxisObjPtr axis = NULL;
      tokens = msStringSplit(value, ',', &num);
      for (j = 0; j < num; ++j) {
        char axisName[500];
        double scale;

        if (msWCSParseScaleString20(tokens[j], axisName, sizeof(axisName),
                                    &scale) != MS_SUCCESS) {
          msFreeCharArray(tokens, num);
          return MS_FAILURE;
        }

        if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
          if (NULL == (axis = msWCSCreateAxisObj20())) {
            msFreeCharArray(tokens, num);
            return MS_FAILURE;
          }
          axis->name = msStrdup(axisName);
          msWCSInsertAxisObj20(params, axis);
        }

        /* check if the size of the axis is already set */
        if (axis->scale != MS_WCS20_UNBOUNDED) {
          msFreeCharArray(tokens, num);
          msSetError(MS_WCSERR, "The scale of the axis is already set.",
                     "msWCSParseRequest20()");
          return MS_FAILURE;
        }
        axis->scale = scale;
      }
      msFreeCharArray(tokens, num);
    } else if (EQUAL(key, "SCALESIZE")) {
      wcs20AxisObjPtr axis = NULL;
      tokens = msStringSplit(value, ',', &num);
      for (j = 0; j < num; ++j) {
        char axisName[500];
        int size;

        if (msWCSParseSizeString20(tokens[j], axisName, sizeof(axisName),
                                   &size) != MS_SUCCESS) {
          msFreeCharArray(tokens, num);
          return MS_FAILURE;
        }

        if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
          if (NULL == (axis = msWCSCreateAxisObj20())) {
            msFreeCharArray(tokens, num);
            return MS_FAILURE;
          }
          axis->name = msStrdup(axisName);
          msWCSInsertAxisObj20(params, axis);
        }

        /* check if the size of the axis is already set */
        if (axis->size != 0) {
          msFreeCharArray(tokens, num);
          msSetError(MS_WCSERR, "The size of the axis is already set.",
                     "msWCSParseRequest20()");
          return MS_FAILURE;
        }
        axis->size = size;
      }
      msFreeCharArray(tokens, num);
    } else if (EQUAL(key, "SCALEEXTENT")) {
      wcs20AxisObjPtr axis = NULL;
      /* No real support for scaleextent, we just interpret it as SCALESIZE */
      tokens = msStringSplit(value, ',', &num);
      for (j = 0; j < num; ++j) {
        char axisName[500];
        int min, max;

        if (msWCSParseScaleExtentString20(tokens[j], axisName, sizeof(axisName),
                                          &min, &max) != MS_SUCCESS) {
          msFreeCharArray(tokens, num);
          return MS_FAILURE;
        }

        if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
          if (NULL == (axis = msWCSCreateAxisObj20())) {
            msFreeCharArray(tokens, num);
            return MS_FAILURE;
          }
          axis->name = msStrdup(axisName);
          msWCSInsertAxisObj20(params, axis);
        }

        /* check if the size of the axis is already set */
        if (axis->size != 0) {
          msFreeCharArray(tokens, num);
          msSetError(MS_WCSERR, "The size of the axis is already set.",
                     "msWCSParseRequest20()");
          return MS_FAILURE;
        }
        axis->size = max - min;
      }
      msFreeCharArray(tokens, num);
      /* We explicitly don't test for strict equality as the parameter name is
       * supposed to be unique */
    } else if (EQUALN(key, "SIZE", 4)) {
      /* Deprecated scaling */
      wcs20AxisObjPtr axis = NULL;
      char axisName[500];
      int size = 0;

      if (msWCSParseSizeString20(value, axisName, sizeof(axisName), &size) ==
          MS_FAILURE) {
        return MS_FAILURE;
      }

      if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
        if (NULL == (axis = msWCSCreateAxisObj20())) {
          return MS_FAILURE;
        }
        axis->name = msStrdup(axisName);
        msWCSInsertAxisObj20(params, axis);
      }

      /* check if the size of the axis is already set */
      if (axis->size != 0) {
        msSetError(MS_WCSERR, "The size of the axis is already set.",
                   "msWCSParseRequest20()");
        return MS_FAILURE;
      }
      axis->size = size;
      /* We explicitly don't test for strict equality as the parameter name is
       * supposed to be unique */
    } else if (EQUALN(key, "RESOLUTION", 10)) {
      wcs20AxisObjPtr axis = NULL;
      char axisName[500];
      double resolution = 0;

      if (msWCSParseResolutionString20(value, axisName, sizeof(axisName),
                                       &resolution) == MS_FAILURE) {
        return MS_FAILURE;
      }

      /* check if axis object already exists, otherwise create a new one */
      if (NULL == (axis = msWCSFindAxis20(params, axisName))) {
        if (NULL == (axis = msWCSCreateAxisObj20())) {
          return MS_FAILURE;
        }
        axis->name = msStrdup(axisName);
        msWCSInsertAxisObj20(params, axis);
      }

      /* check if the resolution of the axis is already set */
      if (axis->resolution != MS_WCS20_UNBOUNDED) {
        msSetError(MS_WCSERR, "The resolution of the axis is already set.",
                   "msWCSParseRequest20()");
        return MS_FAILURE;
      }
      axis->resolution = resolution;
      /* We explicitly don't test for strict equality as the parameter name is
       * supposed to be unique */
    } else if (EQUALN(key, "SUBSET", 6)) {
      wcs20AxisObjPtr axis = NULL;
      wcs20SubsetObjPtr subset = msWCSCreateSubsetObj20();
      if (NULL == subset) {
        return MS_FAILURE;
      }
      if (msWCSParseSubsetKVPString20(subset, value) == MS_FAILURE) {
        msWCSFreeSubsetObj20(subset);
        msWCSException(map, "InvalidSubsetting", "subset",
                       ows_request->version);
        return MS_DONE;
      }

      if (NULL == (axis = msWCSFindAxis20(params, subset->axis))) {
        if (NULL == (axis = msWCSCreateAxisObj20())) {
          return MS_FAILURE;
        }
        axis->name = msStrdup(subset->axis);
        msWCSInsertAxisObj20(params, axis);
      }

      if (NULL != axis->subset) {
        msSetError(MS_WCSERR, "The axis '%s' is already subsetted.",
                   "msWCSParseRequest20()", axis->name);
        msWCSFreeSubsetObj20(subset);
        msWCSException(map, "InvalidAxisLabel", "subset", ows_request->version);
        return MS_DONE;
      }
      axis->subset = subset;
    } else if (EQUAL(key, "RANGESUBSET")) {
      tokens = msStringSplit(value, ',', &num);
      for (j = 0; j < num; ++j) {
        params->range_subset = CSLAddString(params->range_subset, tokens[j]);
      }
      msFreeCharArray(tokens, num);
    } else if (EQUALN(key, "GEOTIFF:", 8)) {
      params->format_options =
          CSLAddNameValue(params->format_options, key, value);
    }
    /* Ignore all other parameters here */
  }

  return MS_SUCCESS;
}

#if defined(USE_LIBXML2)

/************************************************************************/
/*                   msWCSValidateAndFindSubsets20()                    */
/*                                                                      */
/*      Iterates over every axis in the parameters and checks if the    */
/*      axis name is in any string list. If found, but there already is */
/*      a sample found for this axis, an Error is returned. Also if no  */
/*      axis is found for a given axis, an error is returned.           */
/************************************************************************/
static int msWCSValidateAndFindAxes20(wcs20ParamsObjPtr params,
                                      wcs20AxisObjPtr outAxes[]) {
  static const int numAxis = 2;
  const char *const validXAxisNames[] = {
      "x",         "xaxis", "x-axis",   "x_axis",   "long", "long_axis",
      "long-axis", "lon",   "lon_axis", "lon-axis", NULL};
  const char *const validYAxisNames[] = {
      "y", "yaxis", "y-axis", "y_axis", "lat", "lat_axis", "lat-axis", NULL};
  const char *const *const validAxisNames[2] = {validXAxisNames,
                                                validYAxisNames};
  int iParamAxis, iAcceptedAxis, iName, i;

  for (i = 0; i < numAxis; ++i) {
    outAxes[i] = NULL;
  }

  /* iterate over all subsets */
  for (iParamAxis = 0; iParamAxis < params->numaxes; ++iParamAxis) {
    int found = 0;

    /* iterate over all given axes */
    for (iAcceptedAxis = 0; iAcceptedAxis < numAxis; ++iAcceptedAxis) {
      /* iterate over all possible names for the current axis */
      for (iName = 0; validAxisNames[iAcceptedAxis][iName] != NULL; ++iName) {
        /* compare axis name with current possible name */
        if (EQUAL(params->axes[iParamAxis]->name,
                  validAxisNames[iAcceptedAxis][iName])) {
          /* if there is already a sample for the axis, throw error */
          if (outAxes[iAcceptedAxis] != NULL) {
            msSetError(MS_WCSERR,
                       "The axis with the name '%s' corresponds "
                       "to the same axis as the subset with the name '%s'.",
                       "msWCSValidateAndFindAxes20()",
                       outAxes[iAcceptedAxis]->name,
                       params->axes[iParamAxis]->name);
            return MS_FAILURE;
          }

          /* if match is found, save it */
          outAxes[iAcceptedAxis] = params->axes[iParamAxis];
          found = 1;
          break;
        }
      }
      if (found) {
        break;
      }
    }

    /* no valid representation for current subset found */
    /* exit and throw error                             */
    if (found == 0) {
      msSetError(MS_WCSERR, "Invalid subset axis '%s'.",
                 "msWCSValidateAndFindAxes20()",
                 params->axes[iParamAxis]->name);
      return MS_FAILURE;
    }
  }
  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSPrepareNamespaces20()                         */
/*                                                                      */
/*      Inserts namespace definitions into the root node of a DOM       */
/*      structure.                                                      */
/************************************************************************/

static void msWCSPrepareNamespaces20(xmlDocPtr pDoc, xmlNodePtr psRootNode,
                                     mapObj *map, int addInspire) {
  xmlNsPtr psXsiNs;
  char *schemaLocation = NULL;
  char *xsi_schemaLocation = NULL;

  xmlSetNs(psRootNode,
           xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_WCS_20_NAMESPACE_URI,
                    BAD_CAST MS_OWSCOMMON_WCS_NAMESPACE_PREFIX));

  xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_OWS_20_NAMESPACE_URI,
           BAD_CAST MS_OWSCOMMON_OWS_NAMESPACE_PREFIX);
  xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_OGC_NAMESPACE_URI,
           BAD_CAST MS_OWSCOMMON_OGC_NAMESPACE_PREFIX);
  xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_W3C_XSI_NAMESPACE_URI,
           BAD_CAST MS_OWSCOMMON_W3C_XSI_NAMESPACE_PREFIX);
  xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_W3C_XLINK_NAMESPACE_URI,
           BAD_CAST MS_OWSCOMMON_W3C_XLINK_NAMESPACE_PREFIX);
  xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_WCS_20_NAMESPACE_URI,
           BAD_CAST MS_OWSCOMMON_WCS_NAMESPACE_PREFIX);
  xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_GML_32_NAMESPACE_URI,
           BAD_CAST MS_OWSCOMMON_GML_NAMESPACE_PREFIX);
  xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_GMLCOV_10_NAMESPACE_URI,
           BAD_CAST MS_OWSCOMMON_GMLCOV_NAMESPACE_PREFIX);
  xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_SWE_20_NAMESPACE_URI,
           BAD_CAST MS_OWSCOMMON_SWE_NAMESPACE_PREFIX);

  if (addInspire) {
    xmlNewNs(psRootNode, BAD_CAST MS_INSPIRE_COMMON_NAMESPACE_URI,
             BAD_CAST MS_INSPIRE_COMMON_NAMESPACE_PREFIX);
    xmlNewNs(psRootNode, BAD_CAST MS_INSPIRE_DLS_NAMESPACE_URI,
             BAD_CAST MS_INSPIRE_DLS_NAMESPACE_PREFIX);
  }

  psXsiNs = xmlSearchNs(pDoc, psRootNode,
                        BAD_CAST MS_OWSCOMMON_W3C_XSI_NAMESPACE_PREFIX);

  schemaLocation = msEncodeHTMLEntities(msOWSGetSchemasLocation(map));
  xsi_schemaLocation = msStrdup(MS_OWSCOMMON_WCS_20_NAMESPACE_URI);
  xsi_schemaLocation = msStringConcatenate(xsi_schemaLocation, " ");
  xsi_schemaLocation = msStringConcatenate(xsi_schemaLocation, schemaLocation);
  xsi_schemaLocation = msStringConcatenate(
      xsi_schemaLocation, MS_OWSCOMMON_WCS_20_SCHEMAS_LOCATION);
  xsi_schemaLocation = msStringConcatenate(xsi_schemaLocation, " ");

  if (addInspire) {
    xsi_schemaLocation = msStringConcatenate(xsi_schemaLocation,
                                             MS_INSPIRE_DLS_NAMESPACE_URI " ");
    xsi_schemaLocation = msStringConcatenate(
        xsi_schemaLocation, msOWSGetInspireSchemasLocation(map));
    xsi_schemaLocation =
        msStringConcatenate(xsi_schemaLocation, MS_INSPIRE_DLS_SCHEMA_LOCATION);
  }

  xmlNewNsProp(psRootNode, psXsiNs, BAD_CAST "schemaLocation",
               BAD_CAST xsi_schemaLocation);

  msFree(schemaLocation);
  msFree(xsi_schemaLocation);
}

/************************************************************************/
/*                   msWCSGetFormatList20()                             */
/*                                                                      */
/*      Copied from mapwcs.c.                                           */
/************************************************************************/

static char *msWCSGetFormatsList20(mapObj *map, layerObj *layer) {
  char *format_list = msStrdup("");
  char **tokens = NULL, **formats = NULL;
  int i, numtokens = 0, numformats;
  char *value;

  /* -------------------------------------------------------------------- */
  /*      Parse from layer metadata.                                      */
  /* -------------------------------------------------------------------- */
  if (layer != NULL &&
      (value = msOWSGetEncodeMetadata(&(layer->metadata), "CO", "formats",
                                      NULL)) != NULL) {
    tokens = msStringSplit(value, ' ', &numtokens);
    msFree(value);
  }

  /* -------------------------------------------------------------------- */
  /*      Parse from map.web metadata.                                    */
  /* -------------------------------------------------------------------- */
  else if ((value = msOWSGetEncodeMetadata(&(map->web.metadata), "CO",
                                           "formats", NULL)) != NULL) {
    tokens = msStringSplit(value, ' ', &numtokens);
    msFree(value);
  }

  /* -------------------------------------------------------------------- */
  /*      Or generate from all configured raster output formats that      */
  /*      look plausible.                                                 */
  /* -------------------------------------------------------------------- */
  else {
    tokens = (char **)msSmallCalloc(map->numoutputformats, sizeof(char *));
    for (i = 0; i < map->numoutputformats; i++) {
      switch (map->outputformatlist[i]->renderer) {
        /* seemingly normal raster format */
      case MS_RENDER_WITH_AGG:
      case MS_RENDER_WITH_RAWDATA:
        tokens[numtokens++] = msStrdup(map->outputformatlist[i]->name);
        break;
        /* rest of formats aren't really WCS compatible */
      default:
        break;
      }
    }
  }

  /* -------------------------------------------------------------------- */
  /*      Convert outputFormatObj names into mime types and remove        */
  /*      duplicates.                                                     */
  /* -------------------------------------------------------------------- */
  numformats = 0;
  formats = (char **)msSmallCalloc(sizeof(char *), numtokens);

  for (i = 0; i < numtokens; i++) {
    int format_i, j;
    const char *mimetype;

    for (format_i = 0; format_i < map->numoutputformats; format_i++) {
      if (EQUAL(map->outputformatlist[format_i]->name, tokens[i]))
        break;
    }

    if (format_i == map->numoutputformats) {
      msDebug("Failed to find outputformat info on format '%s', ignore.\n",
              tokens[i]);
      continue;
    }

    mimetype = map->outputformatlist[format_i]->mimetype;
    if (mimetype == NULL || strlen(mimetype) == 0) {
      msDebug("No mimetime for format '%s', ignoring.\n", tokens[i]);
      continue;
    }

    for (j = 0; j < numformats; j++) {
      if (EQUAL(mimetype, formats[j]))
        break;
    }

    if (j < numformats) {
      msDebug("Format '%s' ignored since mimetype '%s' duplicates another "
              "outputFormatObj.\n",
              tokens[i], mimetype);
      continue;
    }

    formats[numformats++] = msStrdup(mimetype);
  }

  msFreeCharArray(tokens, numtokens);

  /* -------------------------------------------------------------------- */
  /*      Turn mimetype list into comma delimited form for easy use       */
  /*      with xml functions.                                             */
  /* -------------------------------------------------------------------- */
  for (i = 0; i < numformats; i++) {
    if (i > 0) {
      format_list = msStringConcatenate(format_list, (char *)",");
    }
    format_list = msStringConcatenate(format_list, formats[i]);
  }
  msFreeCharArray(formats, numformats);

  return format_list;
}

/************************************************************************/
/*                   msWCSSwapAxes20                                    */
/*                                                                      */
/*      Helper function to determine if a SRS mandates swapped axes.    */
/************************************************************************/

static int msWCSSwapAxes20(char *srs_uri) {
  int srid = 0;

  /* get SRID from the srs uri */
  if (srs_uri != NULL && strlen(srs_uri) > 0) {
    if (sscanf(srs_uri, "http://www.opengis.net/def/crs/EPSG/0/%d", &srid) !=
        EOF)
      ;
    else if (sscanf(srs_uri, "http://www.opengis.net/def/crs/%d", &srid) != EOF)
      ;
    else
      srid = 0;
  }
  if (srid == 0)
    return MS_FALSE;

  return msIsAxisInverted(srid);
}

/************************************************************************/
/*                   msWCSCommon20_CreateBoundedBy()                    */
/*                                                                      */
/*      Inserts the BoundedBy section into an existing DOM structure.   */
/************************************************************************/

static void msWCSCommon20_CreateBoundedBy(wcs20coverageMetadataObjPtr cm,
                                          xmlNsPtr psGmlNs, xmlNodePtr psRoot,
                                          projectionObj *projection,
                                          int swapAxes) {
  xmlNodePtr psBoundedBy, psEnvelope;
  char lowerCorner[100], upperCorner[100], axisLabels[100], uomLabels[100];

  psBoundedBy = xmlNewChild(psRoot, psGmlNs, BAD_CAST "boundedBy", NULL);
  {
    psEnvelope = xmlNewChild(psBoundedBy, psGmlNs, BAD_CAST "Envelope", NULL);
    {
      xmlNewProp(psEnvelope, BAD_CAST "srsName", BAD_CAST cm->srs_uri);

      if (projection->proj != NULL && msProjIsGeographicCRS(projection)) {
        if (swapAxes == MS_FALSE) {
          strlcpy(axisLabels, "long lat", sizeof(axisLabels));
        } else {
          strlcpy(axisLabels, "lat long", sizeof(axisLabels));
        }
        strlcpy(uomLabels, "deg deg", sizeof(uomLabels));
      } else {
        if (swapAxes == MS_FALSE) {
          strlcpy(axisLabels, "x y", sizeof(axisLabels));
        } else {
          strlcpy(axisLabels, "y x", sizeof(axisLabels));
        }
        strlcpy(uomLabels, "m m", sizeof(uomLabels));
      }
      xmlNewProp(psEnvelope, BAD_CAST "axisLabels", BAD_CAST axisLabels);
      xmlNewProp(psEnvelope, BAD_CAST "uomLabels", BAD_CAST uomLabels);
      xmlNewProp(psEnvelope, BAD_CAST "srsDimension", BAD_CAST "2");

      if (swapAxes == MS_FALSE) {
        snprintf(lowerCorner, sizeof(lowerCorner), "%.15g %.15g",
                 cm->extent.minx, cm->extent.miny);
        snprintf(upperCorner, sizeof(upperCorner), "%.15g %.15g",
                 cm->extent.maxx, cm->extent.maxy);
      } else {
        snprintf(lowerCorner, sizeof(lowerCorner), "%.15g %.15g",
                 cm->extent.miny, cm->extent.minx);
        snprintf(upperCorner, sizeof(upperCorner), "%.15g %.15g",
                 cm->extent.maxy, cm->extent.maxx);
      }

      xmlNewChild(psEnvelope, psGmlNs, BAD_CAST "lowerCorner",
                  BAD_CAST lowerCorner);
      xmlNewChild(psEnvelope, psGmlNs, BAD_CAST "upperCorner",
                  BAD_CAST upperCorner);
    }
  }
}

/************************************************************************/
/*                   msWCSCommon20_CreateDomainSet()                    */
/*                                                                      */
/*      Inserts the DomainSet section into an existing DOM structure.   */
/************************************************************************/

static void msWCSCommon20_CreateDomainSet(layerObj *layer,
                                          wcs20coverageMetadataObjPtr cm,
                                          xmlNsPtr psGmlNs, xmlNodePtr psRoot,
                                          projectionObj *projection,
                                          int swapAxes) {
  xmlNodePtr psDomainSet, psGrid, psLimits, psGridEnvelope, psOrigin, psOffsetX,
      psOffsetY;
  char low[100], high[100], id[100], point[100];
  char offsetVector1[100], offsetVector2[100], axisLabels[100];

  psDomainSet = xmlNewChild(psRoot, psGmlNs, BAD_CAST "domainSet", NULL);
  {
    psGrid = xmlNewChild(psDomainSet, psGmlNs, BAD_CAST "RectifiedGrid", NULL);
    {
      double x0 = cm->geotransform[0] + cm->geotransform[1] / 2 +
                  cm->geotransform[2] / 2;
      double y0 = cm->geotransform[3] + cm->geotransform[4] / 2 +
                  cm->geotransform[5] / 2;
      double resx = cm->geotransform[1];
      double resy = cm->geotransform[5];

      xmlNewProp(psGrid, BAD_CAST "dimension", BAD_CAST "2");
      snprintf(id, sizeof(id), "grid_%s", layer->name);
      xmlNewNsProp(psGrid, psGmlNs, BAD_CAST "id", BAD_CAST id);

      psLimits = xmlNewChild(psGrid, psGmlNs, BAD_CAST "limits", NULL);
      {
        psGridEnvelope =
            xmlNewChild(psLimits, psGmlNs, BAD_CAST "GridEnvelope", NULL);
        {
          strlcpy(low, "0 0", sizeof(low));
          snprintf(high, sizeof(high), "%d %d", cm->xsize - 1, cm->ysize - 1);

          xmlNewChild(psGridEnvelope, psGmlNs, BAD_CAST "low", BAD_CAST low);
          xmlNewChild(psGridEnvelope, psGmlNs, BAD_CAST "high", BAD_CAST high);
        }
      }

      if (projection->proj != NULL && msProjIsGeographicCRS(projection)) {
        strlcpy(axisLabels, "long lat", sizeof(axisLabels));
      } else {
        strlcpy(axisLabels, "x y", sizeof(axisLabels));
      }

      xmlNewChild(psGrid, psGmlNs, BAD_CAST "axisLabels", BAD_CAST axisLabels);

      psOrigin = xmlNewChild(psGrid, psGmlNs, BAD_CAST "origin", NULL);
      {
        if (swapAxes == MS_FALSE) {
          snprintf(point, sizeof(point), "%f %f", x0, y0);
        } else {
          snprintf(point, sizeof(point), "%f %f", y0, x0);
        }
        psOrigin = xmlNewChild(psOrigin, psGmlNs, BAD_CAST "Point", NULL);
        snprintf(id, sizeof(id), "grid_origin_%s", layer->name);
        xmlNewNsProp(psOrigin, psGmlNs, BAD_CAST "id", BAD_CAST id);
        xmlNewProp(psOrigin, BAD_CAST "srsName", BAD_CAST cm->srs_uri);

        xmlNewChild(psOrigin, psGmlNs, BAD_CAST "pos", BAD_CAST point);
      }

      if (swapAxes == MS_FALSE) {
        snprintf(offsetVector1, sizeof(offsetVector1), "%f 0", resx);
        snprintf(offsetVector2, sizeof(offsetVector2), "0 %f", resy);
      } else {
        snprintf(offsetVector1, sizeof(offsetVector1), "0 %f", resx);
        snprintf(offsetVector2, sizeof(offsetVector2), "%f 0", resy);
      }
      psOffsetX = xmlNewChild(psGrid, psGmlNs, BAD_CAST "offsetVector",
                              BAD_CAST offsetVector1);
      psOffsetY = xmlNewChild(psGrid, psGmlNs, BAD_CAST "offsetVector",
                              BAD_CAST offsetVector2);

      xmlNewProp(psOffsetX, BAD_CAST "srsName", BAD_CAST cm->srs_uri);
      xmlNewProp(psOffsetY, BAD_CAST "srsName", BAD_CAST cm->srs_uri);
    }
  }
}

/************************************************************************/
/*                   msWCSCommon20_CreateRangeType()                    */
/*                                                                      */
/*      Inserts the RangeType section into an existing DOM structure.   */
/************************************************************************/

static void msWCSCommon20_CreateRangeType(wcs20coverageMetadataObjPtr cm,
                                          char *bands, xmlNsPtr psGmlcovNs,
                                          xmlNsPtr psSweNs, xmlNodePtr psRoot) {
  xmlNodePtr psRangeType, psDataRecord, psField, psQuantity, psUom,
      psConstraint, psAllowedValues = NULL, psNilValues = NULL;
  char **arr = NULL;
  int num = 0;

  if (NULL != bands) {
    arr = msStringSplit(bands, ',', &num);
  }

  psRangeType = xmlNewChild(psRoot, psGmlcovNs, BAD_CAST "rangeType", NULL);
  psDataRecord = xmlNewChild(psRangeType, psSweNs, BAD_CAST "DataRecord", NULL);

  /* iterate over every band */
  for (unsigned i = 0; i < cm->numbands; ++i) {
    /* only add bands that are in the range subset */
    if (NULL != arr && num > 0) {
      int found = MS_FALSE, j;
      for (j = 0; j < num; ++j) {
        int repr = 0;
        if (msStringParseInteger(arr[j], &repr) == MS_SUCCESS &&
            static_cast<unsigned>(repr) == i + 1) {
          found = MS_TRUE;
          break;
        }
      }
      if (found == MS_FALSE) {
        /* ignore this band since it is not in the range subset */
        continue;
      }
    }

    /* add field tag */
    psField = xmlNewChild(psDataRecord, psSweNs, BAD_CAST "field", NULL);

    if (cm->bands[i].name != NULL) {
      xmlNewProp(psField, BAD_CAST "name", BAD_CAST cm->bands[i].name);
    } else {
      xmlNewProp(psField, BAD_CAST "name", BAD_CAST "band");
    }
    /* add Quantity tag */
    psQuantity = xmlNewChild(psField, psSweNs, BAD_CAST "Quantity", NULL);
    if (cm->bands[i].definition != NULL) {
      xmlNewProp(psQuantity, BAD_CAST "definition",
                 BAD_CAST cm->bands[i].definition);
    }
    if (cm->bands[i].description != NULL) {
      xmlNewChild(psQuantity, psSweNs, BAD_CAST "description",
                  BAD_CAST cm->bands[i].description);
    }

    /* if there are given nilvalues -> add them to the first field */
    /* all other fields get a reference to these */
    if (cm->numnilvalues > 0) {
      psNilValues = xmlNewChild(
          xmlNewChild(psQuantity, psSweNs, BAD_CAST "nilValues", NULL), psSweNs,
          BAD_CAST "NilValues", NULL);
      for (unsigned j = 0; j < cm->numnilvalues; ++j) {
        xmlNodePtr psTemp =
            xmlNewChild(psNilValues, psSweNs, BAD_CAST "nilValue",
                        BAD_CAST cm->nilvalues[j]);
        if (j < cm->numnilvalues)
          xmlNewProp(psTemp, BAD_CAST "reason",
                     BAD_CAST cm->nilvalues_reasons[j]);
      }
    } else { /* create an empty nilValues tag */
      xmlNewChild(psQuantity, psSweNs, BAD_CAST "nilValues", NULL);
    }

    psUom = xmlNewChild(psQuantity, psSweNs, BAD_CAST "uom", NULL);
    if (cm->bands[i].uom != NULL) {
      xmlNewProp(psUom, BAD_CAST "code", BAD_CAST cm->bands[i].uom);
    } else {
      xmlNewProp(psUom, BAD_CAST "code", BAD_CAST "W.m-2.Sr-1");
    }

    /* add constraint */
    psConstraint =
        xmlNewChild(psQuantity, psSweNs, BAD_CAST "constraint", NULL);

    {
      char interval[100], significant_figures[100];
      psAllowedValues =
          xmlNewChild(psConstraint, psSweNs, BAD_CAST "AllowedValues", NULL);

      /* Interval */
      snprintf(interval, sizeof(interval), "%.5g %.5g",
               cm->bands[i].interval_min, cm->bands[i].interval_max);
      xmlNewChild(psAllowedValues, psSweNs, BAD_CAST "interval",
                  BAD_CAST interval);

      /* Significant figures */
      snprintf(significant_figures, sizeof(significant_figures), "%d",
               cm->bands[i].significant_figures);
      xmlNewChild(psAllowedValues, psSweNs, BAD_CAST "significantFigures",
                  BAD_CAST significant_figures);
    }
  }
  msFreeCharArray(arr, num);
}

/************************************************************************/
/*                   msWCSWriteDocument20()                             */
/*                                                                      */
/*      Writes out an xml structure to the stream.                      */
/************************************************************************/

static int msWCSWriteDocument20(xmlDocPtr psDoc) {
  xmlChar *buffer = NULL;
  int size = 0;
  msIOContext *context = NULL;
  char *contenttype = NULL;

  if (msIO_needBinaryStdout() == MS_FAILURE) {
    return MS_FAILURE;
  }

  if (EQUAL((char *)xmlDocGetRootElement(psDoc)->name, "RectifiedGridCoverage"))
    contenttype = msStrdup("application/gml+xml");
  else
    contenttype = msStrdup("text/xml");

  msIO_setHeader("Content-Type", "%s; charset=UTF-8", contenttype);
  msIO_sendHeaders();
  msFree(contenttype);

  context = msIO_getHandler(stdout);

  xmlDocDumpFormatMemoryEnc(psDoc, &buffer, &size, "UTF-8", 1);
  msIO_contextWrite(context, buffer, size);
  xmlFree(buffer);

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSWriteFile20()                                 */
/*                                                                      */
/*      Writes an image object to the stream. If multipart is set,      */
/*      then content sections are inserted.                             */
/************************************************************************/

static int msWCSWriteFile20(mapObj *map, imageObj *image,
                            wcs20ParamsObjPtr params, int multipart) {
  int status;
  char *filename = NULL;
  char *base_dir = NULL;
  const char *fo_filename;
  int i;

  fo_filename = msGetOutputFormatOption(image->format, "FILENAME", NULL);

  /* -------------------------------------------------------------------- */
  /*      Fetch the driver we will be using and check if it supports      */
  /*      VSIL IO.                                                        */
  /* -------------------------------------------------------------------- */
  if (EQUALN(image->format->driver, "GDAL/", 5)) {
    GDALDriverH hDriver;
    const char *pszExtension = image->format->extension;

    msAcquireLock(TLOCK_GDAL);
    hDriver = GDALGetDriverByName(image->format->driver + 5);
    if (hDriver == NULL) {
      msReleaseLock(TLOCK_GDAL);
      msSetError(MS_MISCERR, "Failed to find %s driver.", "msWCSWriteFile20()",
                 image->format->driver + 5);
      return msWCSException(map, "NoApplicableCode", "mapserv",
                            params->version);
    }

    if (pszExtension == NULL)
      pszExtension = "img.tmp";

    if (msGDALDriverSupportsVirtualIOOutput(hDriver)) {
      base_dir = msTmpFile(map, map->mappath, "/vsimem/wcsout", NULL);
      if (fo_filename)
        filename = msStrdup(CPLFormFilename(base_dir, fo_filename, NULL));
      else
        filename = msStrdup(CPLFormFilename(base_dir, "out", pszExtension));

      /*            CleanVSIDir( "/vsimem/wcsout" ); */

      msReleaseLock(TLOCK_GDAL);
      status = msSaveImage(map, image, filename);
      if (status != MS_SUCCESS) {
        msSetError(MS_MISCERR, "msSaveImage() failed", "msWCSWriteFile20()");
        return msWCSException20(map, "NoApplicableCode", "mapserv",
                                params->version);
      }
    }
    msReleaseLock(TLOCK_GDAL);
  }

  /* -------------------------------------------------------------------- */
  /*      If we weren't able to write data under /vsimem, then we just    */
  /*      output a single "stock" filename.                               */
  /* -------------------------------------------------------------------- */
  if (filename == NULL) {
    msOutputFormatResolveFromImage(map, image);
    if (multipart) {
      msIO_fprintf(stdout, "\r\n--wcs\r\n");
      msIO_fprintf(stdout,
                   "Content-Type: %s\r\n"
                   "Content-Description: coverage data\r\n"
                   "Content-Transfer-Encoding: binary\r\n",
                   MS_IMAGE_MIME_TYPE(map->outputformat));

      if (fo_filename != NULL)
        msIO_fprintf(stdout,
                     "Content-ID: coverage/%s\r\n"
                     "Content-Disposition: INLINE; filename=%s\r\n\r\n",
                     fo_filename, fo_filename);
      else
        msIO_fprintf(stdout,
                     "Content-ID: coverage/wcs.%s\r\n"
                     "Content-Disposition: INLINE\r\n\r\n",
                     MS_IMAGE_EXTENSION(map->outputformat));
    } else {
      msIO_setHeader("Content-Type", "%s",
                     MS_IMAGE_MIME_TYPE(map->outputformat));
      msIO_setHeader("Content-Description", "coverage data");
      msIO_setHeader("Content-Transfer-Encoding", "binary");

      if (fo_filename != NULL) {
        msIO_setHeader("Content-ID", "coverage/%s", fo_filename);
        msIO_setHeader("Content-Disposition", "INLINE; filename=%s",
                       fo_filename);
      } else {
        msIO_setHeader("Content-ID", "coverage/wcs.%s",
                       MS_IMAGE_EXTENSION(map->outputformat));
        msIO_setHeader("Content-Disposition", "INLINE");
      }
      msIO_sendHeaders();
    }

    status = msSaveImage(map, image, NULL);
    if (status != MS_SUCCESS) {
      msSetError(MS_MISCERR, "msSaveImage() failed", "msWCSWriteFile20()");
      return msWCSException(map, "NoApplicableCode", "mapserv",
                            params->version);
    }
    if (multipart)
      msIO_fprintf(stdout, "\r\n--wcs--\r\n");
    return MS_SUCCESS;
  }

  /* -------------------------------------------------------------------- */
  /*      When potentially listing multiple files, we take great care     */
  /*      to identify the "primary" file and list it first.  In fact      */
  /*      it is the only file listed in the coverages document.           */
  /* -------------------------------------------------------------------- */
  {
    char **all_files = CPLReadDir(base_dir);
    int count = CSLCount(all_files);

    if (msIO_needBinaryStdout() == MS_FAILURE)
      return MS_FAILURE;

    msAcquireLock(TLOCK_GDAL);
    for (i = count - 1; i >= 0; i--) {
      const char *this_file = all_files[i];

      if (EQUAL(this_file, ".") || EQUAL(this_file, "..")) {
        all_files = CSLRemoveStrings(all_files, i, 1, NULL);
        continue;
      }

      if (i > 0 && EQUAL(this_file, CPLGetFilename(filename))) {
        all_files = CSLRemoveStrings(all_files, i, 1, NULL);
        all_files = CSLInsertString(all_files, 0, CPLGetFilename(filename));
        i++;
      }
    }

    /* -------------------------------------------------------------------- */
    /*      Dump all the files in the memory directory as mime sections.    */
    /* -------------------------------------------------------------------- */
    count = CSLCount(all_files);

    if (count > 1 && multipart == MS_FALSE) {
      msDebug("msWCSWriteFile20(): force multipart output without gml summary "
              "because we have multiple files in the result.\n");

      multipart = MS_TRUE;
      msIO_setHeader("Content-Type", "multipart/related; boundary=wcs");
      msIO_sendHeaders();
    }

    for (i = 0; i < count; i++) {
      const char *mimetype = NULL;
      VSILFILE *fp;
      unsigned char block[4000];
      int bytes_read;

      if (i == 0 && !EQUAL(MS_IMAGE_MIME_TYPE(map->outputformat), "unknown"))
        mimetype = MS_IMAGE_MIME_TYPE(map->outputformat);

      if (mimetype == NULL)
        mimetype = "application/octet-stream";
      if (multipart) {
        msIO_fprintf(stdout, "\r\n--wcs\r\n");
        msIO_fprintf(stdout,
                     "Content-Type: %s\r\n"
                     "Content-Description: coverage data\r\n"
                     "Content-Transfer-Encoding: binary\r\n"
                     "Content-ID: coverage/%s\r\n"
                     "Content-Disposition: INLINE; filename=%s\r\n\r\n",
                     mimetype, all_files[i], all_files[i]);
      } else {
        msIO_setHeader("Content-Type", "%s", mimetype);
        msIO_setHeader("Content-Description", "coverage data");
        msIO_setHeader("Content-Transfer-Encoding", "binary");
        msIO_setHeader("Content-ID", "coverage/%s", all_files[i]);
        msIO_setHeader("Content-Disposition", "INLINE; filename=%s",
                       all_files[i]);
        msIO_sendHeaders();
      }

      fp = VSIFOpenL(CPLFormFilename(base_dir, all_files[i], NULL), "rb");
      if (fp == NULL) {
        msReleaseLock(TLOCK_GDAL);
        msSetError(MS_MISCERR, "Failed to open %s for streaming to stdout.",
                   "msWCSWriteFile20()", all_files[i]);
        return MS_FAILURE;
      }

      while ((bytes_read = VSIFReadL(block, 1, sizeof(block), fp)) > 0)
        msIO_fwrite(block, 1, bytes_read, stdout);

      VSIFCloseL(fp);

      VSIUnlink(CPLFormFilename(base_dir, all_files[i], NULL));
    }

    msFree(base_dir);
    msFree(filename);
    CSLDestroy(all_files);
    msReleaseLock(TLOCK_GDAL);
    if (multipart)
      msIO_fprintf(stdout, "\r\n--wcs--\r\n");
    return MS_SUCCESS;
  }

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSGetRangesetAxisMetadata20()                   */
/*                                                                      */
/*      Looks up a layers metadata for a specific axis information.     */
/************************************************************************/

static const char *msWCSLookupRangesetAxisMetadata20(hashTableObj *table,
                                                     const char *axis,
                                                     const char *item) {
  char buf[500];
  const char *value;

  if (table == NULL || axis == NULL || item == NULL) {
    return NULL;
  }

  strlcpy(buf, axis, sizeof(buf));
  strlcat(buf, "_", sizeof(buf));
  strlcat(buf, item, sizeof(buf));
  if ((value = msLookupHashTable(table, buf)) != NULL) {
    return value;
  }
  return msOWSLookupMetadata(table, "CO", buf);
}

/************************************************************************/
/*                   msWCSGetCoverageMetadata20()                       */
/*                                                                      */
/*      Inits a coverageMetadataObj. Uses msWCSGetCoverageMetadata()    */
/*      but exchanges the SRS URN by an URI for compliance with 2.0.    */
/************************************************************************/

static int msWCSGetCoverageMetadata20(layerObj *layer,
                                      wcs20coverageMetadataObj *cm) {
  char *srs_uri = NULL;
  memset(cm, 0, sizeof(wcs20coverageMetadataObj));
  if (msCheckParentPointer(layer->map, "map") == MS_FAILURE)
    return MS_FAILURE;

  msOWSGetEPSGProj(&(layer->projection), &(layer->metadata), "CO", MS_TRUE,
                   &(cm->srs_epsg));
  if (!cm->srs_epsg) {
    msOWSGetEPSGProj(&(layer->map->projection), &(layer->map->web.metadata),
                     "CO", MS_TRUE, &cm->srs_epsg);
    if (!cm->srs_epsg) {
      msSetError(MS_WCSERR,
                 "Unable to determine the SRS for this layer, "
                 "no projection defined and no metadata available.",
                 "msWCSGetCoverageMetadata20()");
      return MS_FAILURE;
    }
  }

  /* -------------------------------------------------------------------- */
  /*      Get the SRS in uri format.                                      */
  /* -------------------------------------------------------------------- */
  if ((srs_uri = msOWSGetProjURI(&(layer->projection), &(layer->metadata), "CO",
                                 MS_TRUE)) == NULL) {
    srs_uri = msOWSGetProjURI(&(layer->map->projection),
                              &(layer->map->web.metadata), "CO", MS_TRUE);
  }

  if (srs_uri != NULL) {
    if (strlen(srs_uri) > sizeof(cm->srs_uri) - 1) {
      msSetError(MS_WCSERR, "SRS URI too long!", "msWCSGetCoverageMetadata()");
      return MS_FAILURE;
    }

    strlcpy(cm->srs_uri, srs_uri, sizeof(cm->srs_uri));
    msFree(srs_uri);
  } else {
    cm->srs_uri[0] = '\0';
  }

  /* setup nilvalues */
  cm->numnilvalues = 0;
  cm->nilvalues = NULL;
  cm->nilvalues_reasons = NULL;
  cm->native_format = NULL;

  /* -------------------------------------------------------------------- */
  /*      If we have "virtual dataset" metadata on the layer, then use    */
  /*      that in preference to inspecting the file(s).                   */
  /*      We require extent and either size or resolution.                */
  /* -------------------------------------------------------------------- */
  if (msOWSLookupMetadata(&(layer->metadata), "CO", "extent") != NULL &&
      (msOWSLookupMetadata(&(layer->metadata), "CO", "resolution") != NULL ||
       msOWSLookupMetadata(&(layer->metadata), "CO", "size") != NULL)) {
    const char *value;

    /* get extent */
    cm->extent.minx = 0.0;
    cm->extent.maxx = 0.0;
    cm->extent.miny = 0.0;
    cm->extent.maxy = 0.0;
    if (msOWSGetLayerExtent(layer->map, layer, "CO", &cm->extent) == MS_FAILURE)
      return MS_FAILURE;

    /* get resolution */
    cm->xresolution = 0.0;
    cm->yresolution = 0.0;
    if ((value = msOWSLookupMetadata(&(layer->metadata), "CO", "resolution")) !=
        NULL) {
      char **tokens;
      int n;

      tokens = msStringSplit(value, ' ', &n);
      if (tokens == NULL || n != 2) {
        msSetError(MS_WCSERR,
                   "Wrong number of arguments for wcs|ows_resolution metadata.",
                   "msWCSGetCoverageMetadata20()");
        msFreeCharArray(tokens, n);
        return MS_FAILURE;
      }
      cm->xresolution = atof(tokens[0]);
      cm->yresolution = atof(tokens[1]);
      msFreeCharArray(tokens, n);
    }

    /* get Size (in pixels and lines) */
    cm->xsize = 0;
    cm->ysize = 0;
    if ((value = msOWSLookupMetadata(&(layer->metadata), "CO", "size")) !=
        NULL) {
      char **tokens;
      int n;

      tokens = msStringSplit(value, ' ', &n);
      if (tokens == NULL || n != 2) {
        msSetError(MS_WCSERR,
                   "Wrong number of arguments for wcs|ows_size metadata.",
                   "msWCSGetCoverageMetadata20()");
        msFreeCharArray(tokens, n);
        return MS_FAILURE;
      }
      cm->xsize = atoi(tokens[0]);
      cm->ysize = atoi(tokens[1]);
      msFreeCharArray(tokens, n);
    }

    /* try to compute raster size */
    if (cm->xsize == 0 && cm->ysize == 0 && cm->xresolution != 0.0 &&
        cm->yresolution != 0.0 && cm->extent.minx != cm->extent.maxx &&
        cm->extent.miny != cm->extent.maxy) {
      cm->xsize =
          (int)((cm->extent.maxx - cm->extent.minx) / cm->xresolution + 0.5);
      cm->ysize = (int)fabs(
          (cm->extent.maxy - cm->extent.miny) / cm->yresolution + 0.5);
    }

    /* try to compute raster resolution */
    if ((cm->xresolution == 0.0 || cm->yresolution == 0.0) && cm->xsize != 0 &&
        cm->ysize != 0) {
      cm->xresolution = (cm->extent.maxx - cm->extent.minx) / cm->xsize;
      cm->yresolution = (cm->extent.maxy - cm->extent.miny) / cm->ysize;
    }

    /* do we have information to do anything */
    if (cm->xresolution == 0.0 || cm->yresolution == 0.0 || cm->xsize == 0 ||
        cm->ysize == 0) {
      msSetError(MS_WCSERR,
                 "Failed to collect extent and resolution for WCS coverage "
                 "from metadata for layer '%s'.  Need value wcs|ows_resolution "
                 "or wcs|ows_size values.",
                 "msWCSGetCoverageMetadata20()", layer->name);
      return MS_FAILURE;
    }

    /* compute geotransform */
    cm->geotransform[0] = cm->extent.minx;
    cm->geotransform[1] = cm->xresolution;
    cm->geotransform[2] = 0.0;
    cm->geotransform[3] = cm->extent.maxy;
    cm->geotransform[4] = 0.0;
    cm->geotransform[5] = -fabs(cm->yresolution);

    /* get bands count, or assume 1 if not found */
    cm->numbands = 1;
    if ((value = msOWSLookupMetadata(&(layer->metadata), "CO", "bandcount")) !=
        NULL) {
      int numbands = 0;
      if (msStringParseInteger(value, &numbands) != MS_SUCCESS) {
        return MS_FAILURE;
      }
      cm->numbands = (size_t)numbands;
    }

    cm->bands = static_cast<wcs20rasterbandMetadataObj *>(
        msSmallCalloc(sizeof(wcs20rasterbandMetadataObj), cm->numbands));

    /* get bands type, or assume float if not found */
    cm->imagemode = MS_IMAGEMODE_FLOAT32;
    if ((value = msOWSLookupMetadata(&(layer->metadata), "CO", "imagemode")) !=
        NULL) {
      if (EQUAL(value, "INT16"))
        cm->imagemode = MS_IMAGEMODE_INT16;
      else if (EQUAL(value, "FLOAT32"))
        cm->imagemode = MS_IMAGEMODE_FLOAT32;
      else if (EQUAL(value, "BYTE"))
        cm->imagemode = MS_IMAGEMODE_BYTE;
      else {
        msSetError(MS_WCSERR,
                   "Content of wcs|ows_imagemode (%s) not recognised.  Should "
                   "be one of BYTE, INT16 or FLOAT32.",
                   "msWCSGetCoverageMetadata20()", value);
        return MS_FAILURE;
      }
    }

    if ((value = msOWSLookupMetadata(&(layer->metadata), "CO",
                                     "native_format")) != NULL) {
      cm->native_format = msStrdup(value);
    }

    /* determine nilvalues and reasons */
    {
      int n_nilvalues = 0, n_nilvalues_reasons = 0;
      char **t_nilvalues = NULL, **t_nilvalues_reasons = NULL;
      if ((value = msOWSLookupMetadata(&(layer->metadata), "CO",
                                       "nilvalues")) != NULL) {
        t_nilvalues = msStringSplit(value, ' ', &n_nilvalues);
      } else if ((value = msOWSLookupMetadata(&(layer->metadata), "CO",
                                              "rangeset_nullvalue")) != NULL) {
        t_nilvalues = msStringSplit(value, ' ', &n_nilvalues);
      }

      if ((value = msOWSLookupMetadata(&(layer->metadata), "CO",
                                       "nilvalues_reasons")) != NULL) {
        t_nilvalues_reasons = msStringSplit(value, ' ', &n_nilvalues_reasons);
      }

      if (n_nilvalues > 0) {
        int i;
        cm->numnilvalues = n_nilvalues;
        cm->nilvalues =
            static_cast<char **>(msSmallCalloc(sizeof(char *), n_nilvalues));
        cm->nilvalues_reasons =
            static_cast<char **>(msSmallCalloc(sizeof(char *), n_nilvalues));
        for (i = 0; i < n_nilvalues; ++i) {
          cm->nilvalues[i] = msStrdup(t_nilvalues[i]);
          if (i < n_nilvalues_reasons) {
            cm->nilvalues_reasons[i] = msStrdup(t_nilvalues_reasons[i]);
          } else {
            cm->nilvalues_reasons[i] = msStrdup("");
          }
        }
      }
      msFreeCharArray(t_nilvalues, n_nilvalues);
      msFreeCharArray(t_nilvalues_reasons, n_nilvalues_reasons);
    }

    {
      int num_band_names = 0, j;
      char **band_names = NULL;

      const char *wcs11_band_names_key = "rangeset_axes";
      const char *wcs20_band_names_key = "band_names";

      const char *wcs11_interval_key = "interval";
      const char *wcs20_interval_key = "interval";
      const char *interval_key = NULL;

      const char *significant_figures_key = "significant_figures";

      const char *wcs11_keys[] = {"label", "semantic", "values_type",
                                  "values_semantic", "description"};
      const char *wcs20_keys[] = {"band_name", "band_interpretation",
                                  "band_uom", "band_definition",
                                  "band_description"};
      const char **keys = NULL;

      char **interval_array;
      int num_interval;

      wcs20rasterbandMetadataObj default_values;

      /* Decide whether WCS1.1 or WCS2.0 keys should be used */
      if ((value = msOWSLookupMetadata(&(layer->metadata), "CO",
                                       wcs20_band_names_key)) != NULL) {
        keys = wcs20_keys;
        interval_key = wcs20_interval_key;
        band_names = msStringSplit(value, ' ', &num_band_names);
      } else if ((value = msOWSLookupMetadata(&(layer->metadata), "CO",
                                              wcs11_band_names_key)) != NULL) {
        keys = wcs11_keys;
        interval_key = wcs11_interval_key;
        /* "bands" has a special processing in WCS 1.0. See */
        /* msWCSSetDefaultBandsRangeSetInfo */
        if (EQUAL(value, "bands")) {
          num_band_names = cm->numbands;
          band_names = (char **)msSmallMalloc(sizeof(char *) * num_band_names);
          for (int i = 0; i < num_band_names; i++) {
            char szName[30];
            snprintf(szName, sizeof(szName), "Band%d", i + 1);
            band_names[i] = msStrdup(szName);
          }
        } else {
          /* WARNING: in WCS 1.x,, "rangeset_axes" has never been intended */
          /* to contain the list of band names... This code should probably */
          /* be removed */
          band_names = msStringSplit(value, ' ', &num_band_names);
        }
      }

      /* return with error when number of bands does not match    */
      /* the number of names                                      */
      if (static_cast<size_t>(num_band_names) != cm->numbands &&
          num_band_names != 0) {
        msFreeCharArray(band_names, num_band_names);
        msSetError(MS_WCSERR,
                   "Wrong number of band names given in layer '%s'. "
                   "Expected %d, got %d.",
                   "msWCSGetCoverageMetadata20()", layer->name,
                   (int)cm->numbands, num_band_names);
        return MS_FAILURE;
      }

      /* look up default metadata values */
      for (j = 1; j < 5; ++j) {
        if (keys != NULL) {
          default_values.values[j] =
              (char *)msOWSLookupMetadata(&(layer->metadata), "CO", keys[j]);
        }
      }

      /* set default values for interval and significant figures */
      switch (cm->imagemode) {
      case MS_IMAGEMODE_BYTE:
        default_values.interval_min = 0.;
        default_values.interval_max = 255.;
        default_values.significant_figures = 3;
        break;
      case MS_IMAGEMODE_INT16:
        default_values.interval_min = 0.;
        default_values.interval_max = (double)USHRT_MAX;
        default_values.significant_figures = 5;
        break;
      case MS_IMAGEMODE_FLOAT32:
        default_values.interval_min = -FLT_MAX;
        default_values.interval_max = FLT_MAX;
        default_values.significant_figures = 12;
        break;
      default:
        // other imagemodes are invalid (see above), just keep the compiler
        // happy
        msFreeCharArray(band_names, num_band_names);
        return MS_FAILURE;
        break;
      }

      /* lookup default interval values */
      if (interval_key != NULL &&
          (value = msOWSLookupMetadata(&(layer->metadata), "CO",
                                       interval_key)) != NULL) {
        interval_array = msStringSplit(value, ' ', &num_interval);

        if (num_interval != 2 ||
            msStringParseDouble(interval_array[0],
                                &(default_values.interval_min)) != MS_SUCCESS ||
            msStringParseDouble(interval_array[1],
                                &(default_values.interval_max)) != MS_SUCCESS) {
          msFreeCharArray(band_names, num_band_names);
          msFreeCharArray(interval_array, num_interval);
          msSetError(MS_WCSERR, "Wrong interval format for default axis.",
                     "msWCSGetCoverageMetadata20()");
          return MS_FAILURE;
        }
        msFreeCharArray(interval_array, num_interval);
      }

      /* lookup default value for significant figures */
      if ((value = msOWSLookupMetadata(&(layer->metadata), "CO",
                                       significant_figures_key)) != NULL) {
        if (msStringParseInteger(
                value, &(default_values.significant_figures)) != MS_SUCCESS) {
          msFreeCharArray(band_names, num_band_names);
          msSetError(MS_WCSERR,
                     "Wrong significant figures format "
                     "for default axis.",
                     "msWCSGetCoverageMetadata20()");
          return MS_FAILURE;
        }
      }

      /* iterate over every band */
      for (unsigned i = 0; i < cm->numbands; ++i) {
        cm->bands[i].name = NULL;

        /* look up band metadata or copy defaults */
        if (num_band_names != 0) {
          cm->bands[i].name = msStrdup(band_names[i]);
          for (j = 1; j < 5; ++j) {
            value = (char *)msWCSLookupRangesetAxisMetadata20(
                &(layer->metadata), cm->bands[i].name, keys[j]);
            if (value != NULL) {
              cm->bands[i].values[j] = msStrdup(value);
            } else if (default_values.values[j] != NULL) {
              cm->bands[i].values[j] = msStrdup(default_values.values[j]);
            }
          }
        }

        /* set up interval */
        value = (char *)msWCSLookupRangesetAxisMetadata20(
            &(layer->metadata), cm->bands[i].name, interval_key);
        if (value != NULL) {
          num_interval = 0;
          interval_array = msStringSplit(value, ' ', &num_interval);

          if (num_interval != 2 ||
              msStringParseDouble(interval_array[0],
                                  &(cm->bands[i].interval_min)) != MS_SUCCESS ||
              msStringParseDouble(interval_array[1],
                                  &(cm->bands[i].interval_max)) != MS_SUCCESS) {
            msFreeCharArray(band_names, num_band_names);
            msFreeCharArray(interval_array, num_interval);
            msSetError(MS_WCSERR, "Wrong interval format for axis %s.",
                       "msWCSGetCoverageMetadata20()", cm->bands[i].name);
            return MS_FAILURE;
          }

          msFreeCharArray(interval_array, num_interval);
        } else {
          cm->bands[i].interval_min = default_values.interval_min;
          cm->bands[i].interval_max = default_values.interval_max;
        }

        /* set up significant figures */
        value = (char *)msWCSLookupRangesetAxisMetadata20(
            &(layer->metadata), cm->bands[i].name, significant_figures_key);
        if (value != NULL) {
          if (msStringParseInteger(
                  value, &(cm->bands[i].significant_figures)) != MS_SUCCESS) {
            msFreeCharArray(band_names, num_band_names);
            msSetError(MS_WCSERR,
                       "Wrong significant figures format "
                       "for axis %s.",
                       "msWCSGetCoverageMetadata20()", cm->bands[i].name);
            return MS_FAILURE;
          }
        } else {
          cm->bands[i].significant_figures = default_values.significant_figures;
        }
      }

      msFreeCharArray(band_names, num_band_names);
    }
  } else if (layer->data == NULL) {
    /* no virtual metadata, not ok unless we're talking 1 image, hopefully we
     * can fix that */
    msSetError(
        MS_WCSERR,
        "RASTER Layer with no DATA statement and no WCS virtual dataset "
        "metadata.  Tileindexed raster layers not supported for WCS without "
        "virtual dataset metadata (cm->extent, wcs_res, wcs_size).",
        "msWCSGetCoverageMetadata20()");
    return MS_FAILURE;
  } else {
    /* work from the file (e.g. DATA) */
    GDALDatasetH hDS;
    GDALDriverH hDriver;
    GDALRasterBandH hBand;
    char szPath[MS_MAXPATHLEN];
    const char *driver_short_name, *driver_long_name;

    msGDALInitialize();

    msTryBuildPath3(szPath, layer->map->mappath, layer->map->shapepath,
                    layer->data);
    msAcquireLock(TLOCK_GDAL);
    {
      char **connectionoptions =
          msGetStringListFromHashTable(&(layer->connectionoptions));
      hDS = GDALOpenEx(szPath, GDAL_OF_RASTER, NULL,
                       (const char *const *)connectionoptions, NULL);
      CSLDestroy(connectionoptions);
    }
    if (hDS == NULL) {
      msReleaseLock(TLOCK_GDAL);
      msSetError(MS_IOERR, "%s", "msWCSGetCoverageMetadata20()",
                 CPLGetLastErrorMsg());
      return MS_FAILURE;
    }

    msGetGDALGeoTransform(hDS, layer->map, layer, cm->geotransform);

    cm->xsize = GDALGetRasterXSize(hDS);
    cm->ysize = GDALGetRasterYSize(hDS);

    cm->extent.minx = cm->geotransform[0];
    cm->extent.maxx = cm->geotransform[0] + cm->geotransform[1] * cm->xsize +
                      cm->geotransform[2] * cm->ysize;
    cm->extent.miny = cm->geotransform[3] + cm->geotransform[4] * cm->xsize +
                      cm->geotransform[5] * cm->ysize;
    cm->extent.maxy = cm->geotransform[3];

    cm->xresolution = cm->geotransform[1];
    cm->yresolution = cm->geotransform[5];

    cm->numbands = GDALGetRasterCount(hDS);

    /* TODO nilvalues? */

    if (cm->numbands == 0) {
      msReleaseLock(TLOCK_GDAL);
      msSetError(MS_WCSERR,
                 "Raster file %s has no raster bands.  This cannot be used in "
                 "a layer.",
                 "msWCSGetCoverageMetadata20()", layer->data);
      return MS_FAILURE;
    }

    hBand = GDALGetRasterBand(hDS, 1);
    switch (GDALGetRasterDataType(hBand)) {
    case GDT_Byte:
      cm->imagemode = MS_IMAGEMODE_BYTE;
      break;
    case GDT_Int16:
      cm->imagemode = MS_IMAGEMODE_INT16;
      break;
    default:
      cm->imagemode = MS_IMAGEMODE_FLOAT32;
      break;
    }

    cm->bands = static_cast<wcs20rasterbandMetadataObj *>(
        msSmallCalloc(sizeof(wcs20rasterbandMetadataObj), cm->numbands));

    /* get as much band metadata as possible */
    for (unsigned i = 0; i < cm->numbands; ++i) {
      char bandname[32];
      GDALColorInterp colorInterp;
      hBand = GDALGetRasterBand(hDS, i + 1);
      snprintf(bandname, sizeof(bandname), "band%d", i + 1);
      cm->bands[i].name = msStrdup(bandname);
      colorInterp = GDALGetRasterColorInterpretation(hBand);
      cm->bands[i].interpretation =
          msStrdup(GDALGetColorInterpretationName(colorInterp));
      cm->bands[i].uom = msStrdup(GDALGetRasterUnitType(hBand));
      if (strlen(cm->bands[i].uom) == 0) {
        msFree(cm->bands[i].uom);
        cm->bands[i].uom = NULL;
      }

      /* set up interval and significant figures */
      switch (cm->imagemode) {
      case MS_IMAGEMODE_BYTE:
      case MS_IMAGEMODE_PC256:
        cm->bands[i].interval_min = 0.;
        cm->bands[i].interval_max = 255.;
        cm->bands[i].significant_figures = 3;
        break;
      case MS_IMAGEMODE_INT16:
        cm->bands[i].interval_min = 0.;
        cm->bands[i].interval_max = (double)USHRT_MAX;
        cm->bands[i].significant_figures = 5;
        break;
      case MS_IMAGEMODE_FLOAT32:
        cm->bands[i].interval_min = -FLT_MAX;
        cm->bands[i].interval_max = FLT_MAX;
        cm->bands[i].significant_figures = 12;
        break;
      }
    }

    hDriver = GDALGetDatasetDriver(hDS);
    driver_short_name = GDALGetDriverShortName(hDriver);
    driver_long_name = GDALGetDriverLongName(hDriver);
    /* TODO: improve this, exchange strstr() */
    for (int i = 0; i < layer->map->numoutputformats; ++i) {
      if (strstr(layer->map->outputformatlist[i]->driver, driver_short_name) !=
              NULL ||
          strstr(layer->map->outputformatlist[i]->driver, driver_long_name) !=
              NULL) {
        cm->native_format = msStrdup(layer->map->outputformatlist[i]->mimetype);
        break;
      }
    }

    GDALClose(hDS);
    msReleaseLock(TLOCK_GDAL);
  }

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSClearCoverageMetadata20()                     */
/*                                                                      */
/*      Returns all dynamically allocated memory from a coverage meta-  */
/*      data object.                                                    */
/************************************************************************/

static int msWCSClearCoverageMetadata20(wcs20coverageMetadataObj *cm) {
  msFree(cm->native_format);
  for (unsigned i = 0; i < cm->numnilvalues; ++i) {
    msFree(cm->nilvalues[i]);
    msFree(cm->nilvalues_reasons[i]);
  }
  msFree(cm->nilvalues);
  msFree(cm->nilvalues_reasons);

  for (unsigned i = 0; i < cm->numbands; ++i) {
    for (int j = 0; j < 5; ++j) {
      msFree(cm->bands[i].values[j]);
    }
  }
  msFree(cm->bands);
  msFree(cm->srs_epsg);
  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSException20()                                 */
/*                                                                      */
/*      Writes out an OWS compliant exception.                          */
/************************************************************************/

int msWCSException20(mapObj *map, const char *exceptionCode,
                     const char *locator, const char *version) {
  int size = 0;
  const char *status = "400 Bad Request";
  char *errorString = NULL;
  char *schemasLocation = NULL;
  char *xsi_schemaLocation = NULL;
  char version_string[OWS_VERSION_MAXLEN];

  xmlDocPtr psDoc = NULL;
  xmlNodePtr psRootNode = NULL;
  xmlNodePtr psMainNode = NULL;
  xmlNsPtr psNsOws = NULL;
  xmlNsPtr psNsXsi = NULL;
  xmlChar *buffer = NULL;

  errorString = msGetErrorString("\n");
  schemasLocation = msEncodeHTMLEntities(msOWSGetSchemasLocation(map));

  psDoc = xmlNewDoc(BAD_CAST "1.0");

  psRootNode = xmlNewNode(NULL, BAD_CAST "ExceptionReport");
  psNsOws = xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_OWS_20_NAMESPACE_URI,
                     BAD_CAST MS_OWSCOMMON_OWS_NAMESPACE_PREFIX);
  xmlSetNs(psRootNode, psNsOws);

  psNsXsi = xmlNewNs(psRootNode, BAD_CAST MS_OWSCOMMON_W3C_XSI_NAMESPACE_URI,
                     BAD_CAST MS_OWSCOMMON_W3C_XSI_NAMESPACE_PREFIX);

  /* add attributes to root element */
  xmlNewProp(psRootNode, BAD_CAST "version", BAD_CAST version);
  xmlNewProp(psRootNode, BAD_CAST "xml:lang",
             BAD_CAST msOWSGetLanguage(map, "exception"));

  /* get 2 digits version string: '2.0' */
  msOWSGetVersionString(OWS_2_0_0, version_string);
  version_string[3] = '\0';

  xsi_schemaLocation = msStrdup((char *)psNsOws->href);
  xsi_schemaLocation = msStringConcatenate(xsi_schemaLocation, " ");
  xsi_schemaLocation =
      msStringConcatenate(xsi_schemaLocation, (char *)schemasLocation);
  xsi_schemaLocation = msStringConcatenate(xsi_schemaLocation, "/ows/");
  xsi_schemaLocation = msStringConcatenate(xsi_schemaLocation, version_string);
  xsi_schemaLocation =
      msStringConcatenate(xsi_schemaLocation, "/owsExceptionReport.xsd");

  /* add namespace'd attributes to root element */
  xmlNewNsProp(psRootNode, psNsXsi, BAD_CAST "schemaLocation",
               BAD_CAST xsi_schemaLocation);

  /* add child element */
  psMainNode = xmlNewChild(psRootNode, NULL, BAD_CAST "Exception", NULL);

  /* add attributes to child */
  xmlNewProp(psMainNode, BAD_CAST "exceptionCode", BAD_CAST exceptionCode);

  if (locator != NULL) {
    xmlNewProp(psMainNode, BAD_CAST "locator", BAD_CAST locator);
  }

  if (errorString != NULL) {
    char *errorMessage = msEncodeHTMLEntities(errorString);
    xmlNewChild(psMainNode, NULL, BAD_CAST "ExceptionText",
                BAD_CAST errorMessage);
    msFree(errorMessage);
  }

  /*psRootNode = msOWSCommonExceptionReport(psNsOws, OWS_2_0_0,
          schemasLocation, version, msOWSGetLanguage(map, "exception"),
          exceptionCode, locator, errorMessage);*/

  xmlDocSetRootElement(psDoc, psRootNode);

  if (exceptionCode == NULL) {
    /* Do nothing */
  } else if (EQUAL(exceptionCode, "OperationNotSupported") ||
             EQUAL(exceptionCode, "OptionNotSupported")) {
    status = "501 Not Implemented";
  } else if (EQUAL(exceptionCode, "NoSuchCoverage") ||
             EQUAL(exceptionCode, "emptyCoverageIdList") ||
             EQUAL(exceptionCode, "InvalidAxisLabel") ||
             EQUAL(exceptionCode, "InvalidSubsetting")) {
    status = "404 Not Found";
  }

  msIO_setHeader("Status", "%s", status);
  msIO_setHeader("Content-Type", "text/xml; charset=UTF-8");
  msIO_sendHeaders();

  xmlDocDumpFormatMemoryEnc(psDoc, &buffer, &size, "UTF-8", 1);

  msIO_printf("%s", buffer);

  /* free buffer and the document */
  free(errorString);
  free(schemasLocation);
  free(xsi_schemaLocation);
  xmlFree(buffer);
  xmlFreeDoc(psDoc);

  /* clear error since we have already reported it */
  msResetErrorList();
  return MS_FAILURE;
}

#define MAX_MIMES 20

static int msWCSGetCapabilities20_CreateProfiles(
    mapObj *map, xmlNodePtr psServiceIdentification, xmlNsPtr psOwsNs) {
  xmlNodePtr psProfile, psTmpNode;
  char *available_mime_types[MAX_MIMES];
  int i = 0;

  /* even indices are urls, uneven are mime-types, or null*/
  const char *const urls_and_mime_types[] = {
      MS_WCS_20_PROFILE_CORE,
      NULL,
      MS_WCS_20_PROFILE_KVP,
      NULL,
      MS_WCS_20_PROFILE_POST,
      NULL,
      MS_WCS_20_PROFILE_GML,
      NULL,
      MS_WCS_20_PROFILE_GML_MULTIPART,
      NULL,
      MS_WCS_20_PROFILE_GML_SPECIAL,
      NULL,
      MS_WCS_20_PROFILE_GML_GEOTIFF,
      "image/tiff",
      MS_WCS_20_PROFILE_CRS,
      NULL,
      MS_WCS_20_PROFILE_SCALING,
      NULL,
      MS_WCS_20_PROFILE_RANGESUBSET,
      NULL,
      MS_WCS_20_PROFILE_INTERPOLATION,
      NULL,
      NULL,
      NULL /* guardian */
  };

  /* navigate to node where profiles shall be inserted */
  for (psTmpNode = psServiceIdentification->children; psTmpNode->next != NULL;
       psTmpNode = psTmpNode->next) {
    if (EQUAL((char *)psTmpNode->name, "ServiceTypeVersion") &&
        !EQUAL((char *)psTmpNode->next->name, "ServiceTypeVersion"))
      break;
  }

  /* get list of all available mime types */
  msGetOutputFormatMimeList(map, available_mime_types, MAX_MIMES);

  while (urls_and_mime_types[i] != NULL) {
    const char *mime_type = urls_and_mime_types[i + 1];

    /* check if there is a mime type */
    if (mime_type != NULL) {
      /* check if the mime_type is in the list of outputformats */
      if (CSLPartialFindString(available_mime_types, mime_type) == -1)
        continue;
    }

    /* create a new node and attach it in the tree */
    psProfile = xmlNewNode(psOwsNs, BAD_CAST "Profile");
    xmlNodeSetContent(psProfile, BAD_CAST urls_and_mime_types[i]);
    xmlAddNextSibling(psTmpNode, psProfile);

    psTmpNode = psProfile;
    i += 2;
  }
  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSGetCapabilities20_CoverageSummary()           */
/*                                                                      */
/*      Helper function to create a short summary for a specific        */
/*      coverage.                                                       */
/************************************************************************/

static int msWCSGetCapabilities20_CoverageSummary(xmlDocPtr doc,
                                                  xmlNodePtr psContents,
                                                  layerObj *layer) {
  wcs20coverageMetadataObj cm;
  int status;
  xmlNodePtr psCSummary;

  xmlNsPtr psWcsNs =
      xmlSearchNs(doc, xmlDocGetRootElement(doc), BAD_CAST "wcs");

  status = msWCSGetCoverageMetadata20(layer, &cm);
  if (status != MS_SUCCESS)
    return MS_FAILURE;

  psCSummary =
      xmlNewChild(psContents, psWcsNs, BAD_CAST "CoverageSummary", NULL);
  xmlNewChild(psCSummary, psWcsNs, BAD_CAST "CoverageId", BAD_CAST layer->name);
  xmlNewChild(psCSummary, psWcsNs, BAD_CAST "CoverageSubtype",
              BAD_CAST "RectifiedGridCoverage");

  msWCS_11_20_PrintMetadataLinks(layer, doc, psCSummary);

  msWCSClearCoverageMetadata20(&cm);

  return MS_SUCCESS;
}

/************************************************************************/
/*                          msWCSAddInspireDSID20                       */
/************************************************************************/

static void msWCSAddInspireDSID20(mapObj *map, xmlNsPtr psNsInspireDls,
                                  xmlNsPtr psNsInspireCommon,
                                  xmlNodePtr pDlsExtendedCapabilities) {
  const char *dsid_code =
      msOWSLookupMetadata(&(map->web.metadata), "CO", "inspire_dsid_code");
  const char *dsid_ns =
      msOWSLookupMetadata(&(map->web.metadata), "CO", "inspire_dsid_ns");
  if (dsid_code == NULL) {
    xmlAddChild(
        pDlsExtendedCapabilities,
        xmlNewComment(
            BAD_CAST
            "WARNING: Required metadata \"inspire_dsid_code\" missing"));
  } else {
    int ntokensCode = 0, ntokensNS = 0;
    char **tokensCode;
    char **tokensNS = NULL;
    int i;

    tokensCode = msStringSplit(dsid_code, ',', &ntokensCode);
    if (dsid_ns != NULL)
      tokensNS =
          msStringSplitComplex(dsid_ns, ",", &ntokensNS, MS_ALLOWEMPTYTOKENS);
    if (ntokensNS > 0 && ntokensNS != ntokensCode) {
      xmlAddChild(
          pDlsExtendedCapabilities,
          xmlNewComment(
              BAD_CAST
              "WARNING: \"inspire_dsid_code\" and \"inspire_dsid_ns\" have not "
              "the same number of elements. Ignoring inspire_dsid_ns"));
      msFreeCharArray(tokensNS, ntokensNS);
      tokensNS = NULL;
      ntokensNS = 0;
    }
    for (i = 0; i < ntokensCode; i++) {
      xmlNodePtr pSDSI =
          xmlNewNode(psNsInspireDls, BAD_CAST "SpatialDataSetIdentifier");
      xmlAddChild(pDlsExtendedCapabilities, pSDSI);
      xmlNewTextChild(pSDSI, psNsInspireCommon, BAD_CAST "Code",
                      BAD_CAST tokensCode[i]);
      if (ntokensNS > 0 && tokensNS[i][0] != '\0')
        xmlNewTextChild(pSDSI, psNsInspireCommon, BAD_CAST "Namespace",
                        BAD_CAST tokensNS[i]);
    }
    msFreeCharArray(tokensCode, ntokensCode);
    msFreeCharArray(tokensNS, ntokensNS);
  }
}

/************************************************************************/
/*                   msWCSGetCapabilities20()                           */
/*                                                                      */
/*      Performs the GetCapabilities operation. Writes out a xml        */
/*      structure with server specific information and a summary        */
/*      of all available coverages.                                     */
/************************************************************************/

int msWCSGetCapabilities20(mapObj *map, cgiRequestObj *req,
                           wcs20ParamsObjPtr params,
                           owsRequestObj *ows_request) {
  xmlDocPtr psDoc = NULL; /* document pointer */
  xmlNodePtr psRootNode, psOperationsNode, psNode;
  const char *updatesequence = NULL;
  xmlNsPtr psOwsNs = NULL, psXLinkNs = NULL, psWcsNs = NULL, psCrsNs = NULL,
           psIntNs = NULL;
  char *script_url = NULL, *script_url_encoded = NULL, *format_list = NULL;
  int i;
  xmlDocPtr pInspireTmpDoc = NULL;

  const char *inspire_capabilities =
      msOWSLookupMetadata(&(map->web.metadata), "CO", "inspire_capabilities");

  char *validated_language = msOWSLanguageNegotiation(
      map, "CO", params->accept_languages, CSLCount(params->accept_languages));

  /* -------------------------------------------------------------------- */
  /*      Create document.                                                */
  /* -------------------------------------------------------------------- */
  psDoc = xmlNewDoc(BAD_CAST "1.0");

  psRootNode = xmlNewNode(NULL, BAD_CAST "Capabilities");

  xmlDocSetRootElement(psDoc, psRootNode);

  /* -------------------------------------------------------------------- */
  /*      Name spaces                                                     */
  /* -------------------------------------------------------------------- */

  msWCSPrepareNamespaces20(psDoc, psRootNode, map,
                           inspire_capabilities != NULL);

  /* lookup namespaces */
  psOwsNs = xmlSearchNs(psDoc, psRootNode,
                        BAD_CAST MS_OWSCOMMON_OWS_NAMESPACE_PREFIX);
  psWcsNs = xmlSearchNs(psDoc, psRootNode,
                        BAD_CAST MS_OWSCOMMON_WCS_NAMESPACE_PREFIX);
  xmlSearchNs(psDoc, psRootNode, BAD_CAST MS_OWSCOMMON_GML_NAMESPACE_PREFIX);
  psXLinkNs = xmlSearchNs(psDoc, psRootNode, BAD_CAST "xlink");
  psCrsNs = xmlNewNs(psRootNode, BAD_CAST "http://www.opengis.net/wcs/crs/1.0",
                     BAD_CAST "crs");
  psIntNs = xmlNewNs(psRootNode,
                     BAD_CAST "http://www.opengis.net/wcs/interpolation/1.0",
                     BAD_CAST "int");

  xmlSetNs(psRootNode, psWcsNs);

  xmlNewProp(psRootNode, BAD_CAST "version", BAD_CAST params->version);

  updatesequence =
      msOWSLookupMetadata(&(map->web.metadata), "CO", "updatesequence");
  if (params->updatesequence != NULL) {
    i = msOWSNegotiateUpdateSequence(params->updatesequence, updatesequence);
    if (i == 0) { /* current */
      msSetError(
          MS_WCSERR, "UPDATESEQUENCE parameter (%s) is equal to server (%s)",
          "msWCSGetCapabilities20()", params->updatesequence, updatesequence);
      xmlFreeDoc(psDoc);
      return msWCSException(map, "CurrentUpdateSequence", "updatesequence",
                            params->version);
    }
    if (i > 0) { /* invalid */
      msSetError(
          MS_WCSERR, "UPDATESEQUENCE parameter (%s) is higher than server (%s)",
          "msWCSGetCapabilities20()", params->updatesequence, updatesequence);
      xmlFreeDoc(psDoc);
      return msWCSException(map, "InvalidUpdateSequence", "updatesequence",
                            params->version);
    }
  }
  if (updatesequence != NULL) {
    xmlNewProp(psRootNode, BAD_CAST "updateSequence", BAD_CAST updatesequence);
  }

  /* -------------------------------------------------------------------- */
  /*      Service identification.                                         */
  /* -------------------------------------------------------------------- */
  if (MS_WCS_20_CAPABILITIES_INCLUDE_SECTION(params, "ServiceIdentification")) {
    psNode = xmlAddChild(
        psRootNode, msOWSCommonServiceIdentification(psOwsNs, map, "OGC WCS",
                                                     "2.0.1,1.1.1,1.0.0", "CO",
                                                     validated_language));
    msWCSGetCapabilities20_CreateProfiles(map, psNode, psOwsNs);
  }

  /* Service Provider */
  if (MS_WCS_20_CAPABILITIES_INCLUDE_SECTION(params, "ServiceProvider")) {
    xmlAddChild(psRootNode,
                msOWSCommonServiceProvider(psOwsNs, psXLinkNs, map, "CO",
                                           validated_language));
  }

  /* -------------------------------------------------------------------- */
  /*      Operations metadata.                                            */
  /* -------------------------------------------------------------------- */
  if (MS_WCS_20_CAPABILITIES_INCLUDE_SECTION(params, "OperationsMetadata")) {
    if ((script_url = msOWSGetOnlineResource(map, "CO", "onlineresource",
                                             req)) == NULL ||
        (script_url_encoded = msEncodeHTMLEntities(script_url)) == NULL) {
      xmlFreeDoc(psDoc);
      msSetError(MS_WCSERR, "Server URL not found", "msWCSGetCapabilities20()");
      return msWCSException(map, "NoApplicableCode", "mapserv",
                            params->version);
    }
    free(script_url);

    psOperationsNode =
        xmlAddChild(psRootNode, msOWSCommonOperationsMetadata(psOwsNs));

    /* -------------------------------------------------------------------- */
    /*      GetCapabilities - add Sections and AcceptVersions?              */
    /* -------------------------------------------------------------------- */
    psNode = msOWSCommonOperationsMetadataOperation(
        psOwsNs, psXLinkNs, "GetCapabilities", OWS_METHOD_GETPOST,
        script_url_encoded);

    xmlAddChild(psNode->last->last->last,
                msOWSCommonOperationsMetadataDomainType(
                    OWS_2_0_0, psOwsNs, "Constraint", "PostEncoding", "XML"));
    xmlAddChild(psOperationsNode, psNode);

    /* -------------------------------------------------------------------- */
    /*      DescribeCoverage                                                */
    /* -------------------------------------------------------------------- */
    if (msOWSRequestIsEnabled(map, NULL, "C", "DescribeCoverage", MS_FALSE)) {
      psNode = msOWSCommonOperationsMetadataOperation(
          psOwsNs, psXLinkNs, "DescribeCoverage", OWS_METHOD_GETPOST,
          script_url_encoded);
      xmlAddChild(psNode->last->last->last,
                  msOWSCommonOperationsMetadataDomainType(
                      OWS_2_0_0, psOwsNs, "Constraint", "PostEncoding", "XML"));
      xmlAddChild(psOperationsNode, psNode);
    }

    /* -------------------------------------------------------------------- */
    /*      GetCoverage                                                     */
    /* -------------------------------------------------------------------- */
    if (msOWSRequestIsEnabled(map, NULL, "C", "GetCoverage", MS_FALSE)) {
      psNode = msOWSCommonOperationsMetadataOperation(
          psOwsNs, psXLinkNs, "GetCoverage", OWS_METHOD_GETPOST,
          script_url_encoded);

      xmlAddChild(psNode->last->last->last,
                  msOWSCommonOperationsMetadataDomainType(
                      OWS_2_0_0, psOwsNs, "Constraint", "PostEncoding", "XML"));
      xmlAddChild(psOperationsNode, psNode);
    }
    msFree(script_url_encoded);

    /* -------------------------------------------------------------------- */
    /*      Extended Capabilities for inspire                               */
    /* -------------------------------------------------------------------- */

    if (inspire_capabilities) {
      msIOContext *old_context;
      msIOContext *new_context;
      msIOBuffer *buffer;

      xmlNodePtr pRoot;
      xmlNodePtr pOWSExtendedCapabilities;
      xmlNodePtr pDlsExtendedCapabilities;
      xmlNodePtr pChild;

      xmlNsPtr psInspireCommonNs = xmlSearchNs(
          psDoc, psRootNode, BAD_CAST MS_INSPIRE_COMMON_NAMESPACE_PREFIX);
      xmlNsPtr psInspireDlsNs = xmlSearchNs(
          psDoc, psRootNode, BAD_CAST MS_INSPIRE_DLS_NAMESPACE_PREFIX);

      old_context = msIO_pushStdoutToBufferAndGetOldContext();
      msOWSPrintInspireCommonExtendedCapabilities(
          stdout, map, "CO", OWS_WARN, "foo",
          "xmlns:" MS_INSPIRE_COMMON_NAMESPACE_PREFIX
          "=\"" MS_INSPIRE_COMMON_NAMESPACE_URI "\" "
          "xmlns:" MS_INSPIRE_DLS_NAMESPACE_PREFIX
          "=\"" MS_INSPIRE_DLS_NAMESPACE_URI "\" "
          "xmlns:xsi=\"" MS_OWSCOMMON_W3C_XSI_NAMESPACE_URI "\"",
          validated_language, OWS_WCS);

      new_context = msIO_getHandler(stdout);
      buffer = (msIOBuffer *)new_context->cbData;

      /* Remove spaces between > and < to get properly indented result */
      msXMLStripIndentation((char *)buffer->data);

      pInspireTmpDoc = xmlParseDoc((const xmlChar *)buffer->data);
      pRoot = xmlDocGetRootElement(pInspireTmpDoc);
      xmlReconciliateNs(psDoc, pRoot);

      pOWSExtendedCapabilities =
          xmlNewNode(psOwsNs, BAD_CAST "ExtendedCapabilities");
      xmlAddChild(psOperationsNode, pOWSExtendedCapabilities);

      pDlsExtendedCapabilities =
          xmlNewNode(psInspireDlsNs, BAD_CAST "ExtendedCapabilities");
      xmlAddChild(pOWSExtendedCapabilities, pDlsExtendedCapabilities);

      pChild = pRoot->children;
      while (pChild != NULL) {
        xmlNodePtr pNext = pChild->next;
        xmlUnlinkNode(pChild);
        xmlAddChild(pDlsExtendedCapabilities, pChild);
        pChild = pNext;
      }

      msWCSAddInspireDSID20(map, psInspireDlsNs, psInspireCommonNs,
                            pDlsExtendedCapabilities);

      msIO_restoreOldStdoutContext(old_context);
    }
  }

  /* -------------------------------------------------------------------- */
  /*      Service metadata.                                               */
  /* -------------------------------------------------------------------- */

  if (MS_WCS_20_CAPABILITIES_INCLUDE_SECTION(params, "ServiceMetadata")) {
    xmlNodePtr psExtensionNode = NULL;
    psNode = xmlNewChild(psRootNode, psWcsNs, BAD_CAST "ServiceMetadata", NULL);

    /* Apply default formats */
    msApplyDefaultOutputFormats(map);

    /* Add formats list */
    format_list = msWCSGetFormatsList20(map, NULL);
    msLibXml2GenerateList(psNode, psWcsNs, "formatSupported", format_list, ',');
    msFree(format_list);

    psExtensionNode = xmlNewChild(psNode, psWcsNs, BAD_CAST "Extension", NULL);
    /* interpolations supported */
    {
      /* TODO: use URIs/URLs once they are specified */
      const char *const interpolation_methods[] = {"NEAREST", "AVERAGE",
                                                   "BILINEAR"};
      int i;
      xmlNodePtr imNode = xmlNewChild(psExtensionNode, psIntNs,
                                      BAD_CAST "InterpolationMetadata", NULL);

      for (i = 0; i < 3; ++i) {
        xmlNewChild(imNode, psIntNs, BAD_CAST "InterpolationSupported",
                    BAD_CAST interpolation_methods[i]);
      }
    }

    /* Report the supported CRSs */
    {
      char *crs_list = NULL;
      xmlNodePtr crsMetadataNode =
          xmlNewChild(psExtensionNode, psCrsNs, BAD_CAST "CrsMetadata", NULL);

      if ((crs_list = msOWSGetProjURI(&(map->projection), &(map->web.metadata),
                                      "CO", MS_FALSE)) != NULL) {
        msLibXml2GenerateList(crsMetadataNode, psCrsNs, "crsSupported",
                              crs_list, ' ');
        msFree(crs_list);
      } else {
        /* could not determine list of CRSs */
      }
    }
  }

  /* -------------------------------------------------------------------- */
  /*      Contents section.                                               */
  /* -------------------------------------------------------------------- */
  if (MS_WCS_20_CAPABILITIES_INCLUDE_SECTION(params, "Contents")) {
    psNode = xmlNewChild(psRootNode, psWcsNs, BAD_CAST "Contents", NULL);

    if (ows_request->numlayers == 0) {
      xmlAddChild(psNode, xmlNewComment(BAD_CAST(
                              "WARNING: No WCS layers are enabled. "
                              "Check wcs/ows_enable_request settings.")));
    } else {
      for (i = 0; i < map->numlayers; ++i) {
        layerObj *layer = map->layers[i];
        int status;

        if (!msWCSIsLayerSupported(layer))
          continue;

        if (!msIntegerInArray(layer->index, ows_request->enabled_layers,
                              ows_request->numlayers))
          continue;

        status = msWCSGetCapabilities20_CoverageSummary(psDoc, psNode, layer);
        if (status != MS_SUCCESS) {
          msFree(validated_language);
          xmlFreeDoc(psDoc);
          xmlCleanupParser();
          return msWCSException(map, "Internal", "mapserv", params->version);
        }
      }
    }
  }
  /* -------------------------------------------------------------------- */
  /*      Write out the document and clean up.                            */
  /* -------------------------------------------------------------------- */
  msWCSWriteDocument20(psDoc);
  msFree(validated_language);
  xmlFreeDoc(psDoc);
  if (pInspireTmpDoc)
    xmlFreeDoc(pInspireTmpDoc);
  xmlCleanupParser();
  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSDescribeCoverage20_CoverageDescription()      */
/*                                                                      */
/*      Creates a description of a specific coverage with the three     */
/*      major sections: BoundedBy, DomainSet and RangeType.             */
/************************************************************************/

static int msWCSDescribeCoverage20_CoverageDescription(layerObj *layer,
                                                       xmlDocPtr psDoc,
                                                       xmlNodePtr psRootNode) {
  int status, swapAxes;
  wcs20coverageMetadataObj cm;
  xmlNodePtr psCD;
  xmlNsPtr psWcsNs, psGmlNs, psGmlcovNs, psSweNs, psXLinkNs;
  psWcsNs = psGmlNs = psGmlcovNs = psSweNs = psXLinkNs = NULL;

  psWcsNs = xmlSearchNs(psDoc, xmlDocGetRootElement(psDoc),
                        BAD_CAST MS_OWSCOMMON_WCS_NAMESPACE_PREFIX);
  psGmlNs = xmlSearchNs(psDoc, xmlDocGetRootElement(psDoc),
                        BAD_CAST MS_OWSCOMMON_GML_NAMESPACE_PREFIX);
  psGmlcovNs = xmlSearchNs(psDoc, xmlDocGetRootElement(psDoc),
                           BAD_CAST MS_OWSCOMMON_GMLCOV_NAMESPACE_PREFIX);
  psSweNs = xmlSearchNs(psDoc, xmlDocGetRootElement(psDoc),
                        BAD_CAST MS_OWSCOMMON_SWE_NAMESPACE_PREFIX);
  psXLinkNs = xmlSearchNs(psDoc, xmlDocGetRootElement(psDoc),
                          BAD_CAST MS_OWSCOMMON_W3C_XLINK_NAMESPACE_PREFIX);

  /* -------------------------------------------------------------------- */
  /*      Verify layer is processable.                                    */
  /* -------------------------------------------------------------------- */
  if (msCheckParentPointer(layer->map, "map") == MS_FAILURE)
    return MS_FAILURE;

  if (!msWCSIsLayerSupported(layer))
    return MS_SUCCESS;

  /* -------------------------------------------------------------------- */
  /*      Setup coverage metadata.                                        */
  /* -------------------------------------------------------------------- */
  status = msWCSGetCoverageMetadata20(layer, &cm);
  if (status != MS_SUCCESS)
    return status;

  swapAxes = msWCSSwapAxes20(cm.srs_uri);

  /* fill in bands rangeset info, if required. */
  /* msWCSSetDefaultBandsRangeSetInfo( NULL, &cm, layer ); */

  /* -------------------------------------------------------------------- */
  /*      Create CoverageDescription node.                                */
  /* -------------------------------------------------------------------- */
  psCD = xmlNewChild(psRootNode, psWcsNs, BAD_CAST "CoverageDescription", NULL);
  xmlNewNsProp(psCD, psGmlNs, BAD_CAST "id", BAD_CAST layer->name);

  /* -------------------------------------------------------------------- */
  /*      gml:boundedBy                                                   */
  /* -------------------------------------------------------------------- */
  msWCSCommon20_CreateBoundedBy(&cm, psGmlNs, psCD, &(layer->projection),
                                swapAxes);

  xmlNewChild(psCD, psWcsNs, BAD_CAST "CoverageId", BAD_CAST layer->name);

  /* -------------------------------------------------------------------- */
  /*      gml:domainSet                                                   */
  /* -------------------------------------------------------------------- */
  msWCSCommon20_CreateDomainSet(layer, &cm, psGmlNs, psCD, &(layer->projection),
                                swapAxes);

  /* -------------------------------------------------------------------- */
  /*      gmlcov:rangeType                                                */
  /* -------------------------------------------------------------------- */
  msWCSCommon20_CreateRangeType(&cm, NULL, psGmlcovNs, psSweNs, psCD);

  /* -------------------------------------------------------------------- */
  /*      wcs:ServiceParameters                                           */
  /* -------------------------------------------------------------------- */
  {
    xmlNodePtr psSP;

    psSP = xmlNewChild(psCD, psWcsNs, BAD_CAST "ServiceParameters", NULL);
    xmlNewChild(psSP, psWcsNs, BAD_CAST "CoverageSubtype",
                BAD_CAST "RectifiedGridCoverage");

    /* -------------------------------------------------------------------- */
    /*      SupportedCRS                                                    */
    /* -------------------------------------------------------------------- */
    /* for now, WCS 2.0 does not allow per coverage CRS definitions */
    /*{
      xmlNodePtr psSupportedCrss;
      char *owned_value;

      psSupportedCrss = xmlNewChild(psSP, psWcsNs,
                                    BAD_CAST "SupportedCRSs", NULL);

      if ((owned_value = msOWSGetProjURI(&(layer->projection),
                                         &(layer->metadata), "CO", MS_FALSE)) !=
    NULL) { } else if ((owned_value = msOWSGetProjURI(&(layer->map->projection),
                                              &(layer->map->web.metadata), "CO",
    MS_FALSE)) != NULL) { } else { msDebug("missing required information, no
    SRSs defined.\n");
      }

      if (owned_value != NULL && strlen(owned_value) > 0) {
        msLibXml2GenerateList(psSupportedCrss, psWcsNs,
                              "SupportedCRS", owned_value, ' ');
      }

      xmlNewChild(psSupportedCrss, psWcsNs,
                  BAD_CAST "NativeCRS", BAD_CAST cm.srs_uri);

      msFree(owned_value);
    }*/

    /* -------------------------------------------------------------------- */
    /*      SupportedFormats                                                */
    /* -------------------------------------------------------------------- */
    /* for now, WCS 2.0 does not allow per coverage format definitions */
    /*{
      xmlNodePtr psSupportedFormats;
      char *format_list;

      psSupportedFormats =
      xmlNewChild(psSP, psWcsNs, BAD_CAST "SupportedFormats", NULL);

      format_list = msWCSGetFormatsList20(layer->map, layer);

      if (strlen(format_list) > 0) {
        msLibXml2GenerateList(psSupportedFormats, psWcsNs,
                              "SupportedFormat", format_list, ',');

      msFree(format_list);
    }*/

    /* -------------------------------------------------------------------- */
    /*      nativeFormat                                                    */
    /* -------------------------------------------------------------------- */
    xmlNewChild(psSP, psWcsNs, BAD_CAST "nativeFormat",
                BAD_CAST(cm.native_format ? cm.native_format : ""));

    if (!cm.native_format) {
      msDebug("msWCSDescribeCoverage20_CoverageDescription(): "
              "No native format specified.\n");
    }
  }

  msWCSClearCoverageMetadata20(&cm);

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSDescribeCoverage20()                          */
/*                                                                      */
/*      Implementation of the DescribeCoverage Operation. The result     */
/*      of this operation is a xml document, containing specific        */
/*      information about a coverage identified by an ID. The result    */
/*      is written on the stream.                                       */
/************************************************************************/

int msWCSDescribeCoverage20(mapObj *map, wcs20ParamsObjPtr params,
                            owsRequestObj *ows_request) {
  xmlDocPtr psDoc = NULL; /* document pointer */
  xmlNodePtr psRootNode;
  xmlNsPtr psWcsNs = NULL;
  int i, j;

  /* create DOM document and root node */
  psDoc = xmlNewDoc(BAD_CAST "1.0");
  psRootNode = xmlNewNode(NULL, BAD_CAST "CoverageDescriptions");
  xmlDocSetRootElement(psDoc, psRootNode);

  /* prepare initial namespace definitions */
  msWCSPrepareNamespaces20(psDoc, psRootNode, map, MS_FALSE);

  psWcsNs = xmlSearchNs(psDoc, psRootNode,
                        BAD_CAST MS_OWSCOMMON_WCS_NAMESPACE_PREFIX);
  xmlSetNs(psRootNode, psWcsNs);

  /* check if IDs are given */
  if (params->ids) {
    /* for each given ID in the ID-list */
    for (j = 0; params->ids[j]; j++) {
      i = msGetLayerIndex(map, params->ids[j]);
      if (i == -1 || (!msIntegerInArray(GET_LAYER(map, i)->index,
                                        ows_request->enabled_layers,
                                        ows_request->numlayers))) {
        msSetError(MS_WCSERR, "Unknown coverage: (%s)",
                   "msWCSDescribeCoverage20()", params->ids[j]);
        return msWCSException(map, "NoSuchCoverage", "coverage",
                              params->version);
      }
      /* create coverage description for the specified layer */
      if (msWCSDescribeCoverage20_CoverageDescription(
              (GET_LAYER(map, i)), psDoc, psRootNode) == MS_FAILURE) {
        msSetError(MS_WCSERR, "Error retrieving coverage description.",
                   "msWCSDescribeCoverage20()");
        return msWCSException(map, "MissingParameterValue", "coverage",
                              params->version);
      }
    }
  } else {
    /* Throw error, since IDs are mandatory */
    msSetError(MS_WCSERR, "Missing COVERAGEID parameter.",
               "msWCSDescribeCoverage20()");
    return msWCSException(map, "MissingParameterValue", "coverage",
                          params->version);
  }

  /* write out the DOM document to the stream */
  msWCSWriteDocument20(psDoc);

  /* cleanup */
  xmlFreeDoc(psDoc);
  xmlCleanupParser();

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSGetCoverage_FinalizeParamsObj20()             */
/*                                                                      */
/*      Finalizes a wcs20ParamsObj for a GetCoverage operation. In the  */
/*      process, the params bounding box is adjusted to the subsets,     */
/*      width, height and resolution are determined and the subset crs  */
/*      is found out.                                                   */
/************************************************************************/

static int msWCSGetCoverage20_FinalizeParamsObj(wcs20ParamsObjPtr params,
                                                wcs20AxisObjPtr *axes) {
  char *crs = NULL;
  int have_scale, have_size, have_resolution;
  int have_global_scale = (params->scale != MS_WCS20_UNBOUNDED) ? 1 : 0;

  if (axes[0] != NULL) {
    if (axes[0]->subset != NULL) {
      msDebug("Subset for X-axis found: %s\n", axes[0]->subset->axis);
      if (!axes[0]->subset->min.unbounded)
        params->bbox.minx = axes[0]->subset->min.scalar;
      if (!axes[0]->subset->max.unbounded)
        params->bbox.maxx = axes[0]->subset->max.scalar;
      crs = axes[0]->subset->crs;
    }
    params->width = axes[0]->size;
    params->resolutionX = axes[0]->resolution;
    params->scaleX = axes[0]->scale;

    if (axes[0]->resolutionUOM != NULL) {
      params->resolutionUnits = msStrdup(axes[0]->resolutionUOM);
    }

    have_scale = (params->scaleX != MS_WCS20_UNBOUNDED) ? 1 : 0;
    have_size = (params->width != 0) ? 1 : 0;
    have_resolution = (params->resolutionX != MS_WCS20_UNBOUNDED) ? 1 : 0;

    if ((have_global_scale + have_scale + have_size + have_resolution) > 1) {
      msSetError(
          MS_WCSERR,
          "Axis '%s' defines scale, size and/or resolution multiple times.",
          "msWCSGetCoverage20_FinalizeParamsObj()", axes[0]->name);
      return MS_FAILURE;
    }
  }

  if (axes[1] != NULL) {
    if (axes[1]->subset != NULL) {
      msDebug("Subset for Y-axis found: %s\n", axes[1]->subset->axis);
      if (!axes[1]->subset->min.unbounded)
        params->bbox.miny = axes[1]->subset->min.scalar;
      if (!axes[1]->subset->max.unbounded)
        params->bbox.maxy = axes[1]->subset->max.scalar;
      if (crs != NULL && axes[0] != NULL && axes[0]->subset != NULL) {
        if (!EQUAL(crs, axes[1]->subset->crs)) {
          msSetError(MS_WCSERR, "CRS for axis %s and axis %s are not the same.",
                     "msWCSGetCoverage20_FinalizeParamsObj()", axes[0]->name,
                     axes[1]->name);
          return MS_FAILURE;
        }
      } else {
        crs = axes[1]->subset->crs;
      }
    }
    params->height = axes[1]->size;
    params->resolutionY = axes[1]->resolution;
    params->scaleY = axes[1]->scale;

    if (params->resolutionUnits == NULL && axes[1]->resolutionUOM != NULL) {
      params->resolutionUnits = msStrdup(axes[1]->resolutionUOM);
    } else if (params->resolutionUnits != NULL &&
               axes[1]->resolutionUOM != NULL &&
               !EQUAL(params->resolutionUnits, axes[1]->resolutionUOM)) {
      msSetError(MS_WCSERR,
                 "The units of measure of the resolution for"
                 "axis %s and axis %s are not the same.",
                 "msWCSGetCoverage20_FinalizeParamsObj()", axes[0]->name,
                 axes[1]->name);
      return MS_FAILURE;
    }

    have_scale = (params->scaleY != MS_WCS20_UNBOUNDED) ? 1 : 0;
    have_size = (params->height != 0) ? 1 : 0;
    have_resolution = (params->resolutionY != MS_WCS20_UNBOUNDED) ? 1 : 0;

    if ((have_global_scale + have_scale + have_size + have_resolution) > 1) {
      msSetError(
          MS_WCSERR,
          "Axis '%s' defines scale, size and/or resolution multiple times.",
          "msWCSGetCoverage20_FinalizeParamsObj()", axes[1]->name);
      return MS_FAILURE;
    }
  }

  /* if scale is globally set, use this value instead */
  if (params->scale != MS_WCS20_UNBOUNDED) {
    params->scaleX = params->scaleY = params->scale;
  }

  /* check if projections are equal */
  if (crs != NULL) {
    if (params->subsetcrs && !EQUAL(crs, params->subsetcrs)) {
      /* already set and not equal -> raise exception */
      msSetError(MS_WCSERR,
                 "SubsetCRS does not match the CRS of the axes: "
                 "'%s' != '%s'",
                 "msWCSCreateBoundingBox20()", params->subsetcrs, crs);
      return MS_FAILURE;
    } else if (params->subsetcrs) {
      /* equal, and already set -> do nothing */
    } else {
      /* not yet globally set -> set it */
      params->subsetcrs = msStrdup(crs);
    }
  } else if (!params->subsetcrs) {
    /* default to CRS of image */
    /* leave params->subsetcrs to null */
  }

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSGetCoverage20_GetBands()                      */
/*                                                                      */
/*      Returns a string, containing a comma-separated list of band     */
/*      indices.                                                        */
/************************************************************************/

static int msWCSGetCoverage20_GetBands(layerObj *layer,
                                       wcs20ParamsObjPtr params,
                                       wcs20coverageMetadataObjPtr cm,
                                       char **bandlist) {
  int maxlen, index;
  char *interval_stop;
  char **band_ids = NULL;

  /* if rangesubset parameter is not given, default to all bands */
  if (NULL == params->range_subset) {
    *bandlist = msStrdup("1");
    for (unsigned i = 1; i < cm->numbands; ++i) {
      char strnumber[12];
      snprintf(strnumber, sizeof(strnumber), ",%d", i + 1);
      *bandlist = msStringConcatenate(*bandlist, strnumber);
    }
    return MS_SUCCESS;
  }

  maxlen = cm->numbands * 4 * sizeof(char);
  *bandlist = static_cast<char *>(msSmallCalloc(sizeof(char), maxlen));

  /* Use WCS 2.0 metadata items in priority */
  {
    char *tmp =
        msOWSGetEncodeMetadata(&layer->metadata, "CO", "band_names", NULL);

    if (NULL == tmp) {
      /* Otherwise default to WCS 1.x*/
      tmp =
          msOWSGetEncodeMetadata(&layer->metadata, "CO", "rangeset_axes", NULL);
      /* "bands" has a special processing in WCS 1.0. See */
      /* msWCSSetDefaultBandsRangeSetInfo */
      if (tmp != NULL && EQUAL(tmp, "bands")) {
        int num_band_names = cm->numbands;
        band_ids = (char **)msSmallCalloc(sizeof(char *), (num_band_names + 1));
        for (int i = 0; i < num_band_names; i++) {
          char szName[30];
          snprintf(szName, sizeof(szName), "Band%d", i + 1);
          band_ids[i] = msStrdup(szName);
        }
      }
    }

    if (NULL != tmp && band_ids == NULL) {
      band_ids = CSLTokenizeString2(tmp, " ", 0);
    }
    msFree(tmp);
  }

  /* If we still don't have band names, use the band names from the coverage
   * metadata */
  if (band_ids == NULL) {
    band_ids = (char **)CPLCalloc(sizeof(char *), (cm->numbands + 1));
    for (unsigned i = 0; i < cm->numbands; ++i) {
      band_ids[i] = CPLStrdup(cm->bands[i].name);
    }
  }

  /* Iterate over all supplied range */
  const int count = CSLCount(params->range_subset);
  for (int i = 0; i < count; ++i) {
    /* RangeInterval case: defined as "<start>:<stop>" */
    if ((interval_stop = strchr(params->range_subset[i], ':')) != NULL) {
      int start, stop;
      *interval_stop = '\0';
      ++interval_stop;

      /* check if the string represents an integer */
      if (msStringParseInteger(params->range_subset[i], &start) == MS_SUCCESS) {
      }
      /* get index of start clause, or raise an error if none was found */
      else if ((start = CSLFindString(band_ids, params->range_subset[i])) !=
               -1) {
        start += 1; /* adjust index, since GDAL bands start with index 1 */
      } else {
        msSetError(MS_WCSERR, "'%s' is not a valid band identifier.",
                   "msWCSGetCoverage20_GetBands()", params->range_subset[i]);
        return MS_FAILURE;
      }

      /* check if the string represents an integer */
      if (msStringParseInteger(interval_stop, &stop) == MS_SUCCESS) {
      }
      /* get index of stop clause, or raise an error if none was found */
      else if ((stop = CSLFindString(band_ids, interval_stop)) != -1) {
        stop += 1; /* adjust index, since GDAL bands start with index 1 */
      } else {
        msSetError(MS_WCSERR, "'%s' is not a valid band identifier.",
                   "msWCSGetCoverage20_GetBands()", interval_stop);
        return MS_FAILURE;
      }

      /* Check whether or not start and stop are different and stop is higher
       * than start */
      if (stop <= start) {
        msSetError(MS_WCSERR, "Invalid range interval given.",
                   "msWCSGetCoverage20_GetBands()");
        return MS_FAILURE;
      } else if (start < 1 || static_cast<unsigned>(stop) > cm->numbands) {
        msSetError(MS_WCSERR, "Band interval is out of the valid range: 1-%d",
                   "msWCSGetCoverage20_GetBands()", (int)cm->numbands);
        return MS_FAILURE;
      }

      /* expand the interval to a list of indices and push them on the list */
      for (int j = start; j <= stop; ++j) {
        char *tmp;
        if (i != 0 || j != start) {
          strlcat(*bandlist, ",", maxlen);
        }
        tmp = msIntToString(j);
        strlcat(*bandlist, tmp, maxlen);
        msFree(tmp);
      }
    }
    /* RangeComponent case */
    else {
      if (i != 0) {
        strlcat(*bandlist, ",", maxlen);
      }

      /* check if the string represents an integer */
      if (msStringParseInteger(params->range_subset[i], &index) == MS_SUCCESS) {
        char *tmp;
        if (index < 1 || static_cast<unsigned>(index) > cm->numbands) {
          msSetError(MS_WCSERR, "Band index is out of the valid range: 1-%d",
                     "msWCSGetCoverage20_GetBands()", (int)cm->numbands);
          return MS_FAILURE;
        }

        tmp = msIntToString((int)index);
        strlcat(*bandlist, tmp, maxlen);
        msFree(tmp);
        continue;
      }

      /* check if the string is equal to a band identifier    */
      /* if so, what is the index of the band                 */
      if ((index = CSLFindString(band_ids, params->range_subset[i])) != -1) {
        char *tmp = msIntToString((int)index + 1);
        strlcat(*bandlist, tmp, maxlen);
        msFree(tmp);
      } else {
        msSetError(MS_WCSERR, "'%s' is not a valid band identifier.",
                   "msWCSGetCoverage20_GetBands()", params->range_subset[i]);
        return MS_FAILURE;
      }
    }
  }

  CSLDestroy(band_ids);
  return MS_SUCCESS;
}

/* Fixes CPLParseNameValue to allow ':' characters for namespaced keys */

static const char *fixedCPLParseNameValue(const char *string, char **key) {
  const char *value;
  size_t size;
  value = strchr(string, '=');

  if (value == NULL) {
    *key = NULL;
    return NULL;
  }

  size = value - string;
  *key = static_cast<char *>(msSmallMalloc(size + 1));
  strncpy(*key, string, size + 1);
  (*key)[size] = '\0';
  return value + 1;
}

/************************************************************************/
/*                   msWCSSetFormatParams20()                           */
/*                                                                      */
/*      Parses the given format options and sets the appropriate        */
/*      output format values.                                           */
/************************************************************************/

static int msWCSSetFormatParams20(outputFormatObj *format,
                                  char **format_options) {
  /* currently geotiff only */
  char *format_option;
  int i = 0;
  int is_geotiff = (format->mimetype && EQUAL(format->mimetype, "image/tiff"));

  if (!is_geotiff || !format_options) {
    /* Currently only geotiff available */
    return MS_SUCCESS;
  }

  format_option = format_options[0];
  while (format_option) {
    /* key, value */
    char *key;
    const char *value = fixedCPLParseNameValue(format_option, &key);

    if (!key) {
      msSetError(MS_WCSERR, "Could not deduct the option key.",
                 "msWCSSetFormatParams20()");
      return MS_FAILURE;
    } else if (!value) {
      msSetError(MS_WCSERR, "Missing value for parameter '%s'.",
                 "msWCSSetFormatParams20()", key);
      return MS_FAILURE;
    }

    if (EQUAL(key, "geotiff:compression") && is_geotiff) {
      /*COMPRESS=[JPEG/LZW/PACKBITS/DEFLATE/CCITTRLE/CCITTFAX3/CCITTFAX4/NONE]*/
      if (EQUAL(value, "None")) {
        msSetOutputFormatOption(format, "COMPRESS", "NONE");
      } else if (EQUAL(value, "PackBits")) {
        msSetOutputFormatOption(format, "COMPRESS", "PACKBITS");
      } else if (EQUAL(value, "Deflate")) {
        msSetOutputFormatOption(format, "COMPRESS", "DEFLATE");
      } else if (EQUAL(value, "Huffman")) {
        msSetOutputFormatOption(format, "COMPRESS", "CCITTRLE");
      } else if (EQUAL(value, "LZW")) {
        msSetOutputFormatOption(format, "COMPRESS", "LZW");
      } else if (EQUAL(value, "JPEG")) {
        msSetOutputFormatOption(format, "COMPRESS", "JPEG");
      }
      /* unsupported compression methods: CCITTFAX3/CCITTFAX4 */
      else {
        msSetError(MS_WCSERR, "Compression method '%s' not supported.",
                   "msWCSSetFormatParams20()", value);
        return MS_FAILURE;
      }
    }

    else if (EQUAL(key, "geotiff:jpeg_quality") && is_geotiff) {
      int quality;
      /* JPEG_QUALITY=[1-100] */
      if (msStringParseInteger(value, &quality) != MS_SUCCESS) {
        msSetError(MS_WCSERR, "Could not parse jpeg_quality value.",
                   "msWCSSetFormatParams20()");
        return MS_FAILURE;
      } else if (quality < 1 || quality > 100) {
        msSetError(MS_WCSERR, "Invalid jpeg_quality value '%d'.",
                   "msWCSSetFormatParams20()", quality);
        return MS_FAILURE;
      }
      msSetOutputFormatOption(format, "JPEG_QUALITY", value);
    } else if (EQUAL(key, "geotiff:predictor") && is_geotiff) {
      /* PREDICTOR=[None/Horizontal/FloatingPoint] */
      const char *predictor;
      if (EQUAL(value, "None") || EQUAL(value, "1")) {
        predictor = "1";
      } else if (EQUAL(value, "Horizontal") || EQUAL(value, "2")) {
        predictor = "2";
      } else if (EQUAL(value, "FloatingPoint") || EQUAL(value, "3")) {
        predictor = "3";
      } else {
        msSetError(MS_WCSERR, "Invalid predictor value '%s'.",
                   "msWCSSetFormatParams20()", value);
        return MS_FAILURE;
      }
      msSetOutputFormatOption(format, "PREDICTOR", predictor);
    } else if (EQUAL(key, "geotiff:interleave") && is_geotiff) {
      /* INTERLEAVE=[BAND,PIXEL] */
      if (EQUAL(value, "Band")) {
        msSetOutputFormatOption(format, "INTERLEAVE", "BAND");
      } else if (EQUAL(value, "Pixel")) {
        msSetOutputFormatOption(format, "INTERLEAVE", "PIXEL");
      } else {
        msSetError(MS_WCSERR, "Interleave method '%s' not supported.",
                   "msWCSSetFormatParams20()", value);
        return MS_FAILURE;
      }
    } else if (EQUAL(key, "geotiff:tiling") && is_geotiff) {
      /* TILED=YES */
      if (EQUAL(value, "true")) {
        msSetOutputFormatOption(format, "TILED", "YES");
      } else if (EQUAL(value, "false")) {
        msSetOutputFormatOption(format, "TILED", "NO");
      } else {
        msSetError(MS_WCSERR, "Invalid boolean value '%s'.",
                   "msWCSSetFormatParams20()", value);
        return MS_FAILURE;
      }
    } else if (EQUAL(key, "geotiff:tileheight") && is_geotiff) {
      /* BLOCKXSIZE=n */
      int tileheight;

      if (msStringParseInteger(value, &tileheight) != MS_SUCCESS) {
        msSetError(MS_WCSERR, "Could not parse tileheight value.",
                   "msWCSSetFormatParams20()");
        return MS_FAILURE;
      } else if (tileheight < 1 || tileheight % 16) {
        msSetError(MS_WCSERR,
                   "Invalid tileheight value '%d'. "
                   "Must be greater than 0 and dividable by 16.",
                   "msWCSSetFormatParams20()", tileheight);
        return MS_FAILURE;
      }
      msSetOutputFormatOption(format, "BLOCKXSIZE", value);
    } else if (EQUAL(key, "geotiff:tilewidth") && is_geotiff) {
      /* BLOCKYSIZE=n */
      int tilewidth;

      if (msStringParseInteger(value, &tilewidth) != MS_SUCCESS) {
        msSetError(MS_WCSERR, "Could not parse tilewidth value.",
                   "msWCSSetFormatParams20()");
        return MS_FAILURE;
      } else if (tilewidth < 1 || tilewidth % 16) {
        msSetError(MS_WCSERR,
                   "Invalid tilewidth value '%d'. "
                   "Must be greater than 0 and dividable by 16.",
                   "msWCSSetFormatParams20()", tilewidth);
        return MS_FAILURE;
      }
      msSetOutputFormatOption(format, "BLOCKYSIZE", value);
    } else if (EQUALN(key, "geotiff:", 8)) {
      msSetError(MS_WCSERR, "Unrecognized GeoTIFF parameter '%s'.",
                 "msWCSSetFormatParams20()", key);
      return MS_FAILURE;
    }

    msFree(key);

    /* fetch next option */
    format_option = format_options[++i];
  }

  return MS_SUCCESS;
}

/************************************************************************/
/*                   msWCSGetCoverage20()                               */
/*                                                                      */
/*      Implementation of the GetCoverage Operation. The coverage       */
/*      is either returned as an image or as a multipart xml/image.     */
/*      The result is written on the stream.                            */
/************************************************************************/

int msWCSGetCoverage20(mapObj *map, cgiRequestObj *request,
                       wcs20ParamsObjPtr params, owsRequestObj *ows_request) {
  (void)request;
  layerObj *layer = NULL;
  wcs20coverageMetadataObj cm;
  imageObj *image = NULL;
  outputFormatObj *format = NULL;

  rectObj subsets, bbox;
  projectionObj imageProj;

  int status, i;
  double x_1, x_2, y_1, y_2;
  char *coverageName, *bandlist = NULL, numbands[12];

  int doDrawRasterLayerDraw = MS_TRUE;
  GDALDatasetH hDS = NULL;

  int widthFromComputationInImageCRS = 0;
  int heightFromComputationInImageCRS = 0;

  /* number of coverage ids should be 1 */
  if (params->ids == NULL || params->ids[0] == NULL) {
    msSetError(MS_WCSERR, "Required parameter CoverageID was not supplied.",
               "msWCSGetCoverage20()");
    return msWCSException(map, "MissingParameterValue", "coverage",
                          params->version);
  }
  if (params->ids[1] != NULL) {
    msSetError(MS_WCSERR, "GetCoverage operation supports only one coverage.",
               "msWCSGetCoverage20()");
    return msWCSException(map, "TooManyParameterValues", "coverage",
                          params->version);
  }

  /* find the right layer */
  layer = NULL;
  for (i = 0; i < map->numlayers; i++) {
    coverageName = msOWSGetEncodeMetadata(&(GET_LAYER(map, i)->metadata), "CO",
                                          "name", GET_LAYER(map, i)->name);
    if (EQUAL(coverageName, params->ids[0]) &&
        (msIntegerInArray(GET_LAYER(map, i)->index, ows_request->enabled_layers,
                          ows_request->numlayers))) {
      layer = GET_LAYER(map, i);
      i = map->numlayers; /* to exit loop don't use break, we want to free
                             resources first */
    }
    msFree(coverageName);
  }

  /* throw exception if no Layer was found */
  if (layer == NULL) {
    msSetError(
        MS_WCSERR,
        "COVERAGEID=%s not found, not in supported layer list. A layer might be disabled for \
this request. Check wcs/ows_enable_request settings.",
        "msWCSGetCoverage20()", params->ids[0]);
    return msWCSException(map, "NoSuchCoverage", "coverageid", params->version);
  }
  /* retrieve coverage metadata  */
  status = msWCSGetCoverageMetadata20(layer, &cm);
  if (status != MS_SUCCESS) {
    msWCSClearCoverageMetadata20(&cm);
    return MS_FAILURE;
  }

  /* fill in bands rangeset info, if required.  */
  /* msWCSSetDefaultBandsRangeSetInfo(NULL, &cm, layer ); */

  /* set  resolution, size and maximum extent */
  layer->extent = map->extent = cm.extent;
  map->cellsize = cm.xresolution;
  map->width = cm.xsize;
  map->height = cm.ysize;

  /************************************************************************/
  /*      finalize the params object. determine subset crs and subset     */
  /*      bbox. Also project the image to the subset crs.                 */
  /************************************************************************/

  msInitProjection(&imageProj);
  msProjectionInheritContextFrom(&imageProj, &(layer->projection));
  if (msLoadProjectionString(&imageProj, cm.srs_epsg) == -1) {
    msFreeProjection(&imageProj);
    msWCSClearCoverageMetadata20(&cm);
    msSetError(MS_WCSERR, "Error loading CRS %s.", "msWCSGetCoverage20()",
               cm.srs_epsg);
    return msWCSException(map, "InvalidParameterValue", "projection",
                          params->version);
  }

  /* iterate over all subsets and check if they are valid*/
  for (i = 0; i < params->numaxes; ++i) {
    if (params->axes[i]->subset != NULL) {
      if (params->axes[i]->subset->timeOrScalar == MS_WCS20_TIME_VALUE) {
        msFreeProjection(&imageProj);
        msWCSClearCoverageMetadata20(&cm);
        msSetError(MS_WCSERR, "Time values for subsets are not supported. ",
                   "msWCSGetCoverage20()");
        return msWCSException(map, "InvalidSubsetting", "subset",
                              params->version);
      }
      if (params->axes[i]->subset->operation == MS_WCS20_SLICE) {
        msFreeProjection(&imageProj);
        msWCSClearCoverageMetadata20(&cm);
        msSetError(MS_WCSERR, "Subset operation 'slice' is not supported.",
                   "msWCSGetCoverage20()");
        return msWCSException(map, "InvalidSubsetting", "subset",
                              params->version);
      }
    }
  }

  {
    wcs20AxisObjPtr axes[2];
    if (msWCSValidateAndFindAxes20(params, axes) == MS_FAILURE) {
      msFreeProjection(&imageProj);
      msWCSClearCoverageMetadata20(&cm);
      return msWCSException(map, "InvalidAxisLabel", "subset", params->version);
    }
    if (msWCSGetCoverage20_FinalizeParamsObj(params, axes) == MS_FAILURE) {
      msFreeProjection(&imageProj);
      msWCSClearCoverageMetadata20(&cm);
      return msWCSException(map, "InvalidParameterValue", "extent",
                            params->version);
    }
  }

  subsets = params->bbox;

  /* if no subsetCRS was specified use the coverages CRS
     (Requirement 27 of the WCS 2.0 specification) */
  if (!params->subsetcrs) {
    params->subsetcrs = msStrdup(cm.srs_epsg);
  }

  if (EQUAL(params->subsetcrs, "imageCRS")) {
    /* subsets are in imageCRS; reproject them to real coordinates */
    rectObj orig_bbox = subsets;

    msFreeProjection(&(map->projection));
    map->projection = imageProj;

    if (subsets.minx != -DBL_MAX || subsets.maxx != DBL_MAX) {
      x_1 = cm.geotransform[0] + orig_bbox.minx * cm.geotransform[1] +
            orig_bbox.miny * cm.geotransform[2];
      // below +1 look suspicious to me (E. Rouault) and probably mean that
      // there are still issues with pixel-center vs pixel-corner convention
      // The resolution computed on wcs_20_post_getcov_trim_x_y_both_1px
      // with that is 11, whereas it should be 10
      x_2 = cm.geotransform[0] + (orig_bbox.maxx + 1) * cm.geotransform[1] +
            (orig_bbox.maxy + 1) * cm.geotransform[2];

      subsets.minx = MS_MIN(x_1, x_2);
      subsets.maxx = MS_MAX(x_1, x_2);
    }
    if (subsets.miny != -DBL_MAX || subsets.maxy != DBL_MAX) {
      // below +1 are suspicious. See bove comment in the minx/maxx cases
      y_1 = cm.geotransform[3] + (orig_bbox.maxx + 1) * cm.geotransform[4] +
            (orig_bbox.maxy + 1) * cm.geotransform[5];
      /*subsets.miny -= cm.geotransform[4]/2 + cm.geotransform[5]/2;*/
      y_2 = cm.geotransform[3] + orig_bbox.minx * cm.geotransform[4] +
            orig_bbox.miny * cm.geotransform[5];

      subsets.miny = MS_MIN(y_1, y_2);
      subsets.maxy = MS_MAX(y_1, y_2);
    }
  } else { /* if crs is not the 'imageCRS' */
    projectionObj subsetProj;

    /* if the subsets have a crs given, project the image extent to it */
    msInitProjection(&subsetProj);
    msProjectionInheritContextFrom(&subsetProj, &(layer->projection));
    if (msLoadProjectionString(&subsetProj, params->subsetcrs) != MS_SUCCESS) {
      msFreeProjection(&subsetProj);
      msFreeProjection(&imageProj);
      msWCSClearCoverageMetadata20(&cm);
      msSetError(MS_WCSERR, "Error loading CRS %s.", "msWCSGetCoverage20()",
                 params->subsetcrs);
      return msWCSException(map, "InvalidParameterValue", "projection",
                            params->version);
    }

    if (msProjectionsDiffer(&imageProj, &subsetProj)) {

      /* Reprojection of source raster extent of (-180,-90,180,90) to any */
      /* projected CRS is going to exhibit strong anomalies. So instead */
      /* do the reverse, project the subset extent to the layer CRS, and */
      /* see how much the subset extent takes with respect to the source */
      /* raster extent. This is only used if output width and resolutionX (or */
      /* (height and resolutionY) are unknown. */
      if (((params->width == 0 && params->resolutionX == MS_WCS20_UNBOUNDED) ||
           (params->height == 0 &&
            params->resolutionY == MS_WCS20_UNBOUNDED)) &&
          (msProjIsGeographicCRS(&imageProj) &&
           !msProjIsGeographicCRS(&subsetProj) &&
           fabs(layer->extent.minx - -180.0) < 1e-5 &&
           fabs(layer->extent.miny - -90.0) < 1e-5 &&
           fabs(layer->extent.maxx - 180.0) < 1e-5 &&
           fabs(layer->extent.maxy - 90.0) < 1e-5)) {
        rectObj subsetInImageProj = subsets;
        if (msProjectRect(&subsetProj, &imageProj, &(subsetInImageProj)) ==
            MS_SUCCESS) {
          subsetInImageProj.minx =
              MS_MAX(subsetInImageProj.minx, layer->extent.minx);
          subsetInImageProj.miny =
              MS_MAX(subsetInImageProj.miny, layer->extent.miny);
          subsetInImageProj.maxx =
              MS_MIN(subsetInImageProj.maxx, layer->extent.maxx);
          subsetInImageProj.maxy =
              MS_MIN(subsetInImageProj.maxy, layer->extent.maxy);
          {
            double total = std::abs(layer->extent.maxx - layer->extent.minx);
            double part =
                std::abs(subsetInImageProj.maxx - subsetInImageProj.minx);
            widthFromComputationInImageCRS =
                MS_NINT((part * map->width) / total);
          }
          {
            double total = std::abs(layer->extent.maxy - layer->extent.miny);
            double part =
                std::abs(subsetInImageProj.maxy - subsetInImageProj.miny);
            heightFromComputationInImageCRS =
                MS_NINT((part * map->height) / total);
          }
        }
      }

      msProjectRect(&imageProj, &subsetProj, &(layer->extent));
      if (msProjIsGeographicCRS(&subsetProj)) {
        if (layer->extent.minx > 180 && subsets.maxx <= 180) {
          layer->extent.minx -= 360;
          layer->extent.maxx -= 360;
        } else if (layer->extent.maxx < -180 && subsets.minx >= -180) {
          layer->extent.minx += 360;
          layer->extent.maxx += 360;
        }
      }
      map->extent = layer->extent;
      msFreeProjection(&(map->projection));
      map->projection = subsetProj;
      msFreeProjection(&imageProj);
    } else {
      msFreeProjection(&(map->projection));
      map->projection = imageProj;
      msFreeProjection(&subsetProj);
    }
  }

  /* create boundings of params subsets and image extent */
  if (msRectOverlap(&subsets, &(layer->extent)) == MS_FALSE) {
    /* extent and bbox do not overlap -> exit */
    msWCSClearCoverageMetadata20(&cm);
    msSetError(MS_WCSERR,
               "Image extent does not intersect with desired region.",
               "msWCSGetCoverage20()");
    return msWCSException(map, "ExtentError", "extent", params->version);
  }

  /* write combined bounding box */
  bbox.minx = MS_MAX(subsets.minx, map->extent.minx);
  bbox.miny = MS_MAX(subsets.miny, map->extent.miny);
  bbox.maxx = MS_MIN(subsets.maxx, map->extent.maxx);
  bbox.maxy = MS_MIN(subsets.maxy, map->extent.maxy);

  /* check if we are overspecified  */
  if ((params->width != 0 && params->resolutionX != MS_WCS20_UNBOUNDED) ||
      (params->height != 0 && params->resolutionY != MS_WCS20_UNBOUNDED)) {
    msWCSClearCoverageMetadata20(&cm);
    msSetError(MS_WCSERR,
               "GetCoverage operation supports only one of SIZE or RESOLUTION "
               "per axis.",
               "msWCSGetCoverage20()");
    return msWCSException(map, "TooManyParameterValues", "coverage",
                          params->version);
  }

  /************************************************************************/
  /* check both axes: see if either size or resolution are given (and     */
  /* calculate the other value). If both are not given, calculate them    */
  /* from the bounding box.                                               */
  /************************************************************************/

  /* check x axis */
  if (params->width != 0) {
    /* TODO Unit Of Measure? */
    params->resolutionX = (bbox.maxx - bbox.minx) / params->width;
  } else if (params->resolutionX != MS_WCS20_UNBOUNDED) {
    params->width = MS_NINT((bbox.maxx - bbox.minx) / params->resolutionX);
  } else {
    if (widthFromComputationInImageCRS != 0) {
      params->width = widthFromComputationInImageCRS;
    } else if (std::abs(bbox.maxx - bbox.minx) !=
               std::abs(map->extent.maxx - map->extent.minx)) {
      double total = std::abs(map->extent.maxx - map->extent.minx),
             part = std::abs(bbox.maxx - bbox.minx);
      params->width = MS_NINT((part * map->width) / total);
    } else {
      params->width = map->width;
    }

    params->resolutionX = (bbox.maxx - bbox.minx) / params->width;

    if (params->scaleX != MS_WCS20_UNBOUNDED) {
      params->resolutionX /= params->scaleX;
      params->width = (long)(((double)params->width) * params->scaleX);
    }
  }

  /* check y axis */
  if (params->height != 0) {
    params->resolutionY = (bbox.maxy - bbox.miny) / params->height;
  } else if (params->resolutionY != MS_WCS20_UNBOUNDED) {
    params->height = MS_NINT((bbox.maxy - bbox.miny) / params->resolutionY);
  } else {
    if (heightFromComputationInImageCRS != 0) {
      params->height = heightFromComputationInImageCRS;
    } else if (std::abs(bbox.maxy - bbox.miny) !=
               std::abs(map->extent.maxy - map->extent.miny)) {
      double total = std::abs(map->extent.maxy - map->extent.miny),
             part = std::abs(bbox.maxy - bbox.miny);
      params->height = MS_NINT((part * map->height) / total);
    } else {
      params->height = map->height;
    }

    params->resolutionY = (bbox.maxy - bbox.miny) / params->height;

    if (params->scaleY != MS_WCS20_UNBOUNDED) {
      params->resolutionY /= params->scaleY;
      params->height = (long)(((double)params->height) * params->scaleY);
    }
  }

  /* WCS 2.0 is center of pixel oriented */
  rectObj bboxOriginIsCorner = bbox;
  bbox.minx += params->resolutionX * 0.5;
  bbox.maxx -= params->resolutionX * 0.5;
  bbox.miny += params->resolutionY * 0.5;
  bbox.maxy -= params->resolutionY * 0.5;

  /* if parameter 'outputcrs' is given, project the image to this crs */
  if (params->outputcrs != NULL) {
    projectionObj outputProj;

    msInitProjection(&outputProj);
    msProjectionInheritContextFrom(&outputProj, &(layer->projection));
    if (msLoadProjectionString(&outputProj, params->outputcrs) == -1) {
      msFreeProjection(&outputProj);
      msWCSClearCoverageMetadata20(&cm);
      return msWCSException(map, "InvalidParameterValue", "coverage",
                            params->version);
    }
    if (msProjectionsDiffer(&(map->projection), &outputProj)) {

      msDebug("msWCSGetCoverage20(): projecting to outputcrs %s\n",
              params->outputcrs);

      msProjectRect(&(map->projection), &outputProj, &bbox);

      /* recalculate resolutions, needed if UOM changes (e.g: deg -> m) */
      if (params->width == 1 || params->height == 1) {
        // if one of the dimension is 1, we cannot use the center-of-pixel bbox,
        // but must use the origin-is-corner one.
        msProjectRect(&(map->projection), &outputProj, &bboxOriginIsCorner);
        params->resolutionX =
            (bboxOriginIsCorner.maxx - bboxOriginIsCorner.minx) / params->width;
        params->resolutionY =
            (bboxOriginIsCorner.maxy - bboxOriginIsCorner.miny) /
            params->height;
      } else {
        params->resolutionX = (bbox.maxx - bbox.minx) / (params->width - 1);
        params->resolutionY = (bbox.maxy - bbox.miny) / (params->height - 1);
      }

      msFreeProjection(&(map->projection));
      map->projection = outputProj;
    } else {
      msFreeProjection(&outputProj);
    }
  }

  /* set the bounding box as new map extent */
  map->extent = bbox;
  map->width = params->width;
  map->height = params->height;

  /* Are we exceeding the MAXSIZE limit on result size? */
  if (map->width > map->maxsize || map->height > map->maxsize) {
    msWCSClearCoverageMetadata20(&cm);
    msSetError(MS_WCSERR,
               "Raster size out of range, width and height of "
               "resulting coverage must be no more than MAXSIZE=%d.",
               "msWCSGetCoverage20()", map->maxsize);

    return msWCSException(map, "InvalidParameterValue", "size",
                          params->version);
  }

  /* Mapserver only supports square cells */
  if (params->resolutionX <= params->resolutionY)
    map->cellsize = params->resolutionX;
  else
    map->cellsize = params->resolutionY;

  msDebug("msWCSGetCoverage20(): Set parameters from original"
          "data. Width: %d, height: %d, cellsize: %f, extent: %f,%f,%f,%f\n",
          map->width, map->height, map->cellsize, map->extent.minx,
          map->extent.miny, map->extent.maxx, map->extent.maxy);

  /**
   * Which format to use?
   *
   * 1) format parameter
   * 2) native format (from metadata) or GDAL format of the input dataset
   * 3) exception
   **/

  if (!params->format) {
    if (cm.native_format) {
      params->format = msStrdup(cm.native_format);
    }
  }

  if (!params->format) {
    msSetError(MS_WCSERR,
               "Output format could not be automatically determined. "
               "Use the FORMAT parameter to specify a format.",
               "msWCSGetCoverage20()");
    msWCSClearCoverageMetadata20(&cm);
    return msWCSException(map, "MissingParameterValue", "format",
                          params->version);
  }

  /*    make sure layer is on   */
  layer->status = MS_ON;

  if (params->width == 1 || params->height == 1)
    msMapComputeGeotransformEx(map, params->resolutionX, params->resolutionY);
  else
    msMapComputeGeotransform(map);

  /*    fill in bands rangeset info, if required.  */
  /* msWCSSetDefaultBandsRangeSetInfo(params, &cm, layer); */
  /* msDebug("Bandcount: %d\n", cm.bandcount); */

  msApplyDefaultOutputFormats(map);

  if (msGetOutputFormatIndex(map, params->format) == -1) {
    msSetError(MS_WCSERR, "Unrecognized value '%s' for the FORMAT parameter.",
               "msWCSGetCoverage20()", params->format);
    msWCSClearCoverageMetadata20(&cm);
    return msWCSException(map, "InvalidParameterValue", "format",
                          params->version);
  }

  /* create a temporary outputformat (we likely will need to tweak parts) */
  format = msCloneOutputFormat(msSelectOutputFormat(map, params->format));
  msApplyOutputFormat(&(map->outputformat), format, MS_NOOVERRIDE);

  /* set format specific parameters */
  if (msWCSSetFormatParams20(format, params->format_options) != MS_SUCCESS) {
    msFree(bandlist);
    msWCSClearCoverageMetadata20(&cm);
    return msWCSException(map, "InvalidParameterValue", "format",
                          params->version);
  }

  if (msWCSGetCoverage20_GetBands(layer, params, &cm, &bandlist) !=
      MS_SUCCESS) {
    msFree(bandlist);
    msWCSClearCoverageMetadata20(&cm);
    return msWCSException(map, "InvalidParameterValue", "rangesubset",
                          params->version);
  }
  msLayerSetProcessingKey(layer, "BANDS", bandlist);
  snprintf(numbands, sizeof(numbands), "%d", msCountChars(bandlist, ',') + 1);
  msSetOutputFormatOption(map->outputformat, "BAND_COUNT", numbands);

  msWCSApplyLayerCreationOptions(layer, map->outputformat, bandlist);
  msWCSApplyLayerMetadataItemOptions(layer, map->outputformat, bandlist);

  /* check for the interpolation */
  /* Defaults to NEAREST */
  if (params->interpolation != NULL) {
    if (EQUALN(params->interpolation, "NEAREST", 7)) {
      msLayerSetProcessingKey(layer, "RESAMPLE", "NEAREST");
    } else if (EQUAL(params->interpolation, "BILINEAR")) {
      msLayerSetProcessingKey(layer, "RESAMPLE", "BILINEAR");
    } else if (EQUAL(params->interpolation, "AVERAGE")) {
      msLayerSetProcessingKey(layer, "RESAMPLE", "AVERAGE");
    } else {
      msFree(bandlist);
      msSetError(MS_WCSERR,
                 "'%s' specifies an unsupported interpolation method.",
                 "msWCSGetCoverage20()", params->interpolation);
      msWCSClearCoverageMetadata20(&cm);
      return msWCSException(map, "InvalidParameterValue", "interpolation",
                            params->version);
    }
  } else {
    msLayerSetProcessingKey(layer, "RESAMPLE", "NEAREST");
  }

  /* since the dataset is only used in one layer, set it to be    */
  /* closed after drawing the layer. This normally defaults to    */
  /* DEFER and will produce a memory leak, because the dataset    */
  /* will not be closed.                                          */
  if (msLayerGetProcessingKey(layer, "CLOSE_CONNECTION") == NULL) {
    msLayerSetProcessingKey(layer, "CLOSE_CONNECTION", "NORMAL");
  }

  if (layer->tileindex == NULL && layer->data != NULL &&
      strlen(layer->data) > 0 && layer->connectiontype != MS_KERNELDENSITY) {
    if (msDrawRasterLayerLowCheckIfMustDraw(map, layer)) {
      char *decrypted_path = NULL;
      char szPath[MS_MAXPATHLEN];
      hDS = (GDALDatasetH)msDrawRasterLayerLowOpenDataset(
          map, layer, layer->data, szPath, &decrypted_path);
      msFree(decrypted_path);
      if (hDS) {
        msWCSApplyDatasetMetadataAsCreationOptions(layer, map->outputformat,
                                                   bandlist, hDS);
        msWCSApplySourceDatasetMetadata(layer, map->outputformat, bandlist,
                                        hDS);
      }
    } else {
      doDrawRasterLayerDraw = MS_FALSE;
    }
  }

  /* create the image object  */
  if (MS_RENDERER_PLUGIN(map->outputformat)) {
    image =
        msImageCreate(map->width, map->height, map->outputformat,
                      map->web.imagepath, map->web.imageurl, map->resolution,
                      map->defresolution, &map->imagecolor);
  } else if (MS_RENDERER_RAWDATA(map->outputformat)) {
    image =
        msImageCreate(map->width, map->height, map->outputformat,
                      map->web.imagepath, map->web.imageurl, map->resolution,
                      map->defresolution, &map->imagecolor);
  } else {
    msFree(bandlist);
    msWCSClearCoverageMetadata20(&cm);
    msSetError(MS_WCSERR, "Map outputformat not supported for WCS!",
               "msWCSGetCoverage20()");
    msDrawRasterLayerLowCloseDataset(layer, hDS);
    return msWCSException(map, NULL, NULL, params->version);
  }

  if (image == NULL) {
    msFree(bandlist);
    msWCSClearCoverageMetadata20(&cm);
    msDrawRasterLayerLowCloseDataset(layer, hDS);
    return msWCSException(map, NULL, NULL, params->version);
  }

  if (layer->mask) {
    int maskLayerIdx = msGetLayerIndex(map, layer->mask);
    layerObj *maskLayer;
    outputFormatObj *altFormat;
    if (maskLayerIdx == -1) {
      msSetError(MS_MISCERR, "Layer (%s) references unknown mask layer (%s)",
                 "msDrawLayer()", layer->name, layer->mask);
      msFreeImage(image);
      msFree(bandlist);
      msWCSClearCoverageMetadata20(&cm);
      msDrawRasterLayerLowCloseDataset(layer, hDS);
      return msWCSException(map, NULL, NULL, params->version);
    }
    maskLayer = GET_LAYER(map, maskLayerIdx);
    if (!maskLayer->maskimage) {
      int i, retcode;
      int origstatus, origlabelcache;
      char *origImageType = msStrdup(map->imagetype);
      altFormat = msSelectOutputFormat(map, "png24");
      msInitializeRendererVTable(altFormat);
      /* TODO: check the png24 format hasn't been tampered with, i.e. it's agg
       */
      maskLayer->maskimage = msImageCreate(
          image->width, image->height, altFormat, image->imagepath,
          image->imageurl, map->resolution, map->defresolution, NULL);
      if (!maskLayer->maskimage) {
        msSetError(MS_MISCERR, "Unable to initialize mask image.",
                   "msDrawLayer()");
        msFreeImage(image);
        msFree(bandlist);
        msWCSClearCoverageMetadata20(&cm);
        msDrawRasterLayerLowCloseDataset(layer, hDS);
        return msWCSException(map, NULL, NULL, params->version);
      }

      /*
       * force the masked layer to status on, and turn off the labelcache so
       * that eventual labels are added to the temporary image instead of being
       * added to the labelcache
       */
      origstatus = maskLayer->status;
      origlabelcache = maskLayer->labelcache;
      maskLayer->status = MS_ON;
      maskLayer->labelcache = MS_OFF;

      /* draw the mask layer in the temporary image */
      retcode = msDrawLayer(map, maskLayer, maskLayer->maskimage);
      maskLayer->status = origstatus;
      maskLayer->labelcache = origlabelcache;
      if (retcode != MS_SUCCESS) {
        free(origImageType);
        msFreeImage(image);
        msFree(bandlist);
        msWCSClearCoverageMetadata20(&cm);
        msDrawRasterLayerLowCloseDataset(layer, hDS);
        return msWCSException(map, NULL, NULL, params->version);
      }
      /*
       * hack to work around bug #3834: if we have use an alternate renderer,
       * the symbolset may contain symbols that reference it. We want to remove
       * those references before the altFormat is destroyed to avoid a segfault
       * and/or a leak, and so the the main renderer doesn't pick the cache up
       * thinking it's for him.
       */
      for (i = 0; i < map->symbolset.numsymbols; i++) {
        if (map->symbolset.symbol[i] != NULL) {
          symbolObj *s = map->symbolset.symbol[i];
          if (s->renderer == MS_IMAGE_RENDERER(maskLayer->maskimage)) {
            MS_IMAGE_RENDERER(maskLayer->maskimage)->freeSymbol(s);
            s->renderer = NULL;
          }
        }
      }
      /* set the imagetype from the original outputformat back (it was removed
       * by msSelectOutputFormat() */
      msFree(map->imagetype);
      map->imagetype = origImageType;
    }
  }

  /* Actually produce the "grid". */
  if (MS_RENDERER_RAWDATA(map->outputformat)) {
    if (doDrawRasterLayerDraw) {
      status = msDrawRasterLayerLowWithDataset(map, layer, image, NULL, hDS);
    } else {
      status = MS_SUCCESS;
    }
  } else {
    rasterBufferObj rb;
    status = MS_IMAGE_RENDERER(image)->getRasterBufferHandle(image, &rb);
    if (MS_LIKELY(status == MS_SUCCESS)) {
      if (doDrawRasterLayerDraw) {
        status = msDrawRasterLayerLowWithDataset(map, layer, image, &rb, hDS);
      } else {
        status = MS_SUCCESS;
      }
    }
  }

  msDrawRasterLayerLowCloseDataset(layer, hDS);

  if (status != MS_SUCCESS) {
    msFree(bandlist);
    msFreeImage(image);
    msWCSClearCoverageMetadata20(&cm);
    return msWCSException(map, NULL, NULL, params->version);
  }

  /* GML+Image */
  /* Embed the image into multipart message */
  if (params->multipart == MS_TRUE) {
    xmlDocPtr psDoc = NULL; /* document pointer */
    xmlNodePtr psRootNode, psRangeSet, psFile, psRangeParameters;
    xmlNsPtr psGmlNs = NULL, psGmlcovNs = NULL, psSweNs = NULL,
             psXLinkNs = NULL;
    wcs20coverageMetadataObj tmpCm;
    char *srs_uri, *default_filename;
    const char *filename;
    int swapAxes;

    /* Create Document  */
    psDoc = xmlNewDoc(BAD_CAST "1.0");
    psRootNode = xmlNewNode(
        NULL, BAD_CAST MS_WCS_GML_COVERAGETYPE_RECTIFIED_GRID_COVERAGE);
    xmlDocSetRootElement(psDoc, psRootNode);

    msWCSPrepareNamespaces20(psDoc, psRootNode, map, MS_FALSE);

    psGmlNs = xmlSearchNs(psDoc, psRootNode,
                          BAD_CAST MS_OWSCOMMON_GML_NAMESPACE_PREFIX);
    psGmlcovNs = xmlSearchNs(psDoc, psRootNode,
                             BAD_CAST MS_OWSCOMMON_GMLCOV_NAMESPACE_PREFIX);
    psSweNs = xmlSearchNs(psDoc, psRootNode,
                          BAD_CAST MS_OWSCOMMON_SWE_NAMESPACE_PREFIX);
    xmlSearchNs(psDoc, psRootNode, BAD_CAST MS_OWSCOMMON_WCS_NAMESPACE_PREFIX);
    psXLinkNs = xmlSearchNs(psDoc, psRootNode,
                            BAD_CAST MS_OWSCOMMON_W3C_XLINK_NAMESPACE_PREFIX);

    xmlNewNsProp(psRootNode, psGmlNs, BAD_CAST "id", BAD_CAST layer->name);

    xmlSetNs(psRootNode, psGmlcovNs);

    srs_uri = msOWSGetProjURI(&map->projection, NULL, "CO", 1);

    tmpCm = cm;
    tmpCm.extent = map->extent;
    tmpCm.xsize = map->width;
    tmpCm.ysize = map->height;
    strlcpy(tmpCm.srs_uri, srs_uri, sizeof(tmpCm.srs_uri));

    tmpCm.xresolution = map->gt.geotransform[1];
    tmpCm.yresolution = map->gt.geotransform[5];

    tmpCm.extent.minx =
        MS_MIN(map->gt.geotransform[0],
               map->gt.geotransform[0] + map->width * tmpCm.xresolution);
    tmpCm.extent.miny =
        MS_MIN(map->gt.geotransform[3],
               map->gt.geotransform[3] + map->height * tmpCm.yresolution);
    tmpCm.extent.maxx =
        MS_MAX(map->gt.geotransform[0],
               map->gt.geotransform[0] + map->width * tmpCm.xresolution);
    tmpCm.extent.maxy =
        MS_MAX(map->gt.geotransform[3],
               map->gt.geotransform[3] + map->height * tmpCm.yresolution);

    swapAxes = msWCSSwapAxes20(srs_uri);
    msFree(srs_uri);

    /* Setup layer information  */
    msWCSCommon20_CreateBoundedBy(&tmpCm, psGmlNs, psRootNode,
                                  &(map->projection), swapAxes);
    msWCSCommon20_CreateDomainSet(layer, &tmpCm, psGmlNs, psRootNode,
                                  &(map->projection), swapAxes);

    psRangeSet = xmlNewChild(psRootNode, psGmlNs, BAD_CAST "rangeSet", NULL);
    psFile = xmlNewChild(psRangeSet, psGmlNs, BAD_CAST "File", NULL);

    /* TODO: wait for updated specifications */
    psRangeParameters =
        xmlNewChild(psFile, psGmlNs, BAD_CAST "rangeParameters", NULL);

    default_filename = msStrdup("out.");
    default_filename = msStringConcatenate(default_filename,
                                           MS_IMAGE_EXTENSION(image->format));

    filename =
        msGetOutputFormatOption(image->format, "FILENAME", default_filename);
    std::string file_ref("cid:coverage/");
    file_ref += filename;

    msFree(default_filename);

    std::string role;
    if (EQUAL(MS_IMAGE_MIME_TYPE(map->outputformat), "image/tiff")) {
      role = MS_WCS_20_PROFILE_GML_GEOTIFF;
    } else {
      role = MS_IMAGE_MIME_TYPE(map->outputformat);
    }

    xmlNewNsProp(psRangeParameters, psXLinkNs, BAD_CAST "href",
                 BAD_CAST file_ref.c_str());
    xmlNewNsProp(psRangeParameters, psXLinkNs, BAD_CAST "role",
                 BAD_CAST role.c_str());
    xmlNewNsProp(psRangeParameters, psXLinkNs, BAD_CAST "arcrole",
                 BAD_CAST "fileReference");

    xmlNewChild(psFile, psGmlNs, BAD_CAST "fileReference",
                BAD_CAST file_ref.c_str());
    xmlNewChild(psFile, psGmlNs, BAD_CAST "fileStructure", NULL);
    xmlNewChild(psFile, psGmlNs, BAD_CAST "mimeType",
                BAD_CAST MS_IMAGE_MIME_TYPE(map->outputformat));

    msWCSCommon20_CreateRangeType(&cm, bandlist, psGmlcovNs, psSweNs,
                                  psRootNode);

    msIO_setHeader("Content-Type", "multipart/related; boundary=wcs");
    msIO_sendHeaders();
    msIO_printf("\r\n--wcs\r\n");

    msWCSWriteDocument20(psDoc);
    msWCSWriteFile20(map, image, params, 1);

    xmlFreeDoc(psDoc);
    xmlCleanupParser();
    /* just print out the file without gml */
  } else {
    msWCSWriteFile20(map, image, params, 0);
  }

  msFree(bandlist);
  msWCSClearCoverageMetadata20(&cm);
  msFreeImage(image);
  return MS_SUCCESS;
}

#endif /* defined(USE_LIBXML2) */

#endif /* defined(USE_WCS_SVR) */

Coverage Report

Created: 2026-04-01 06:20