/******************************************************************************
 *
 * Project:  WCS Client Driver
 * Purpose:  Implementation of Dataset class for WCS 2.0.
 * Author:   Ari Jolma <ari dot jolma at gmail dot com>
 *
 ******************************************************************************
 * Copyright (c) 2006, Frank Warmerdam
 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
 * Copyright (c) 2017, Ari Jolma
 * Copyright (c) 2017, Finnish Environment Institute
 *
 * 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 or substantial portions of the 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 "cpl_string.h"
#include "cpl_minixml.h"
#include "cpl_http.h"
#include "cpl_conv.h"
#include "gmlutils.h"
#include "gdal_frmts.h"
#include "gdal_pam.h"
#include "ogr_spatialref.h"
#include "gmlcoverage.h"

#include <algorithm>

#include "wcsdataset.h"
#include "wcsutils.h"

using namespace WCSUtils;

/************************************************************************/
/*                         CoverageSubtype()                            */
/*                                                                      */
/************************************************************************/

static CPLString CoverageSubtype(CPLXMLNode *coverage)
{
    CPLString subtype =
        CPLGetXMLValue(coverage, "ServiceParameters.CoverageSubtype", "");
    size_t pos = subtype.find("Coverage");
    if (pos != std::string::npos)
    {
        subtype.erase(pos);
    }
    return subtype;
}

/************************************************************************/
/*                         GetGridNode()                                */
/*                                                                      */
/************************************************************************/

static CPLXMLNode *GetGridNode(CPLXMLNode *coverage, const CPLString &subtype)
{
    CPLXMLNode *grid = nullptr;
    // Construct the name of the node that we look under domainSet.
    // For now we can handle RectifiedGrid and ReferenceableGridByVectors.
    // Note that if this is called at GetCoverage stage, the grid should not be
    // NULL.
    CPLString path = "domainSet";
    if (subtype == "RectifiedGrid")
    {
        grid = CPLGetXMLNode(coverage, (path + "." + subtype).c_str());
    }
    else if (subtype == "ReferenceableGrid")
    {
        grid = CPLGetXMLNode(coverage,
                             (path + "." + subtype + "ByVectors").c_str());
    }
    if (!grid)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Can't handle coverages of type '%s'.", subtype.c_str());
    }
    return grid;
}

/************************************************************************/
/*                         ParseParameters()                            */
/*                                                                      */
/************************************************************************/

static void ParseParameters(CPLXMLNode *service,
                            std::vector<CPLString> &dimensions,
                            CPLString &range,
                            std::vector<std::vector<CPLString>> &others)
{
    std::vector<CPLString> parameters =
        Split(CPLGetXMLValue(service, "Parameters", ""), "&");
    for (unsigned int i = 0; i < parameters.size(); ++i)
    {
        std::vector<CPLString> kv = Split(parameters[i], "=");
        if (kv.size() < 2)
        {
            continue;
        }
        kv[0].toupper();
        if (kv[0] == "RANGESUBSET")
        {
            range = kv[1];
        }
        else if (kv[0] == "SUBSET")
        {
            dimensions = Split(kv[1], ";");
        }
        else
        {
            std::vector<CPLString> kv2;
            kv2.push_back(kv[0]);
            kv2.push_back(kv[1]);
            others.push_back(kv2);
        }
    }
    // fallback to service values, if any
    if (range == "")
    {
        range = CPLGetXMLValue(service, "RangeSubset", "");
    }
    if (dimensions.size() == 0)
    {
        dimensions = Split(CPLGetXMLValue(service, "Subset", ""), ";");
    }
}

/************************************************************************/
/*                         GetExtent()                                  */
/*                                                                      */
/************************************************************************/

std::vector<double> WCSDataset201::GetExtent(int nXOff, int nYOff, int nXSize,
                                             int nYSize,
                                             CPL_UNUSED int nBufXSize,
                                             CPL_UNUSED int nBufYSize)
{
    std::vector<double> extent;
    // WCS 2.0 extents are the outer edges of outer pixels.
    extent.push_back(adfGeoTransform[0] + (nXOff)*adfGeoTransform[1]);
    extent.push_back(adfGeoTransform[3] +
                     (nYOff + nYSize) * adfGeoTransform[5]);
    extent.push_back(adfGeoTransform[0] +
                     (nXOff + nXSize) * adfGeoTransform[1]);
    extent.push_back(adfGeoTransform[3] + (nYOff)*adfGeoTransform[5]);
    return extent;
}

/************************************************************************/
/*                        GetCoverageRequest()                          */
/*                                                                      */
/************************************************************************/

CPLString WCSDataset201::GetCoverageRequest(bool scaled, int nBufXSize,
                                            int nBufYSize,
                                            const std::vector<double> &extent,
                                            CPL_UNUSED CPLString osBandList)
{
    CPLString request = CPLGetXMLValue(psService, "ServiceURL", "");
    request = CPLURLAddKVP(request, "SERVICE", "WCS");
    request += "&REQUEST=GetCoverage";
    request +=
        "&VERSION=" + CPLString(CPLGetXMLValue(psService, "Version", ""));
    request += "&COVERAGEID=" +
               URLEncode(CPLGetXMLValue(psService, "CoverageName", ""));

    // note: native_crs is not really supported
    if (!native_crs)
    {
        CPLString crs = URLEncode(CPLGetXMLValue(psService, "SRS", ""));
        request += "&OUTPUTCRS=" + crs;
        request += "&SUBSETTINGCRS=" + crs;
    }

    std::vector<CPLString> domain =
        Split(CPLGetXMLValue(psService, "Domain", ""), ",");
    if (domain.size() < 2)
    {
        // eek!
        domain.push_back("E");
        domain.push_back("N");
    }
    const char *x = domain[0].c_str();
    const char *y = domain[1].c_str();
    if (CPLGetXMLBoolean(psService, "SubsetAxisSwap"))
    {
        const char *tmp = x;
        x = y;
        y = tmp;
    }

    std::vector<CPLString> low =
        Split(CPLGetXMLValue(psService, "Low", ""), ",");
    std::vector<CPLString> high =
        Split(CPLGetXMLValue(psService, "High", ""), ",");
    CPLString a = CPLString().Printf("%.17g", extent[0]);
    if (low.size() > 1 && CPLAtof(low[0].c_str()) > extent[0])
    {
        a = low[0];
    }
    CPLString b = CPLString().Printf("%.17g", extent[2]);
    if (high.size() > 1 && CPLAtof(high[0].c_str()) < extent[2])
    {
        b = high[0];
    }
    /*
    CPLString a = CPLString().Printf(
        "%.17g", MAX(adfGeoTransform[0], extent[0]));
    CPLString b = CPLString().Printf(
        "%.17g", MIN(adfGeoTransform[0] + nRasterXSize * adfGeoTransform[1],
    extent[2]));
    */

    // 09-147 KVP Protocol: subset keys must be unique
    // GeoServer: seems to require plain SUBSET for x and y

    request +=
        CPLString().Printf("&SUBSET=%s%%28%s,%s%%29", x, a.c_str(), b.c_str());

    a = CPLString().Printf("%.17g", extent[1]);
    if (low.size() > 1 && CPLAtof(low[1].c_str()) > extent[1])
    {
        a = low[1];
    }
    b = CPLString().Printf("%.17g", extent[3]);
    if (high.size() > 1 && CPLAtof(high[1].c_str()) < extent[3])
    {
        b = high[1];
    }
    /*
    a = CPLString().Printf(
        "%.17g", MAX(adfGeoTransform[3] + nRasterYSize * adfGeoTransform[5],
    extent[1])); b = CPLString().Printf(
        "%.17g", MIN(adfGeoTransform[3], extent[3]));
    */

    request +=
        CPLString().Printf("&SUBSET=%s%%28%s,%s%%29", y, a.c_str(), b.c_str());

    // Dimension and range parameters:
    std::vector<CPLString> dimensions;
    CPLString range;
    std::vector<std::vector<CPLString>> others;
    ParseParameters(psService, dimensions, range, others);

    // set subsets for axis other than x/y
    for (unsigned int i = 0; i < dimensions.size(); ++i)
    {
        size_t pos = dimensions[i].find("(");
        CPLString dim = dimensions[i].substr(0, pos);
        if (IndexOf(dim, domain) != -1)
        {
            continue;
        }
        std::vector<CPLString> params =
            Split(FromParenthesis(dimensions[i]), ",");
        request +=
            "&SUBSET" + CPLString().Printf("%i", i) + "=" + dim + "%28";  // (
        for (unsigned int j = 0; j < params.size(); ++j)
        {
            // todo: %22 (") should be used only for non-numbers
            request += "%22" + params[j] + "%22";
        }
        request += "%29";  // )
    }

    if (scaled)
    {
        CPLString tmp;
        // scaling is expressed in grid axes
        if (CPLGetXMLBoolean(psService, "UseScaleFactor"))
        {
            double fx = fabs((extent[2] - extent[0]) / adfGeoTransform[1] /
                             ((double)nBufXSize + 0.5));
            double fy = fabs((extent[3] - extent[1]) / adfGeoTransform[5] /
                             ((double)nBufYSize + 0.5));
            tmp.Printf("&SCALEFACTOR=%.15g", MIN(fx, fy));
        }
        else
        {
            std::vector<CPLString> grid_axes =
                Split(CPLGetXMLValue(psService, "GridAxes", ""), ",");
            if (grid_axes.size() < 2)
            {
                // eek!
                grid_axes.push_back("E");
                grid_axes.push_back("N");
            }
            tmp.Printf("&SCALESIZE=%s%%28%i%%29,%s%%28%i%%29",
                       grid_axes[0].c_str(), nBufXSize, grid_axes[1].c_str(),
                       nBufYSize);
        }
        request += tmp;
    }

    if (range != "" && range != "*")
    {
        request += "&RANGESUBSET=" + range;
    }

    // other parameters may come from
    // 1) URL (others)
    // 2) Service file
    const char *keys[] = {WCS_URL_PARAMETERS};
    for (unsigned int i = 0; i < CPL_ARRAYSIZE(keys); i++)
    {
        CPLString value;
        int ix = IndexOf(CPLString(keys[i]).toupper(), others);
        if (ix >= 0)
        {
            value = others[ix][1];
        }
        else
        {
            value = CPLGetXMLValue(psService, keys[i], "");
        }
        if (value != "")
        {
            request = CPLURLAddKVP(request, keys[i], value);
        }
    }
    // add extra parameters
    CPLString extra = CPLGetXMLValue(psService, "Parameters", "");
    if (extra != "")
    {
        std::vector<CPLString> pairs = Split(extra, "&");
        for (unsigned int i = 0; i < pairs.size(); ++i)
        {
            std::vector<CPLString> pair = Split(pairs[i], "=");
            request = CPLURLAddKVP(request, pair[0], pair[1]);
        }
    }
    std::vector<CPLString> pairs =
        Split(CPLGetXMLValue(psService, "GetCoverageExtra", ""), "&");
    for (unsigned int i = 0; i < pairs.size(); ++i)
    {
        std::vector<CPLString> pair = Split(pairs[i], "=");
        if (pair.size() > 1)
        {
            request = CPLURLAddKVP(request, pair[0], pair[1]);
        }
    }

    CPLDebug("WCS", "Requesting %s", request.c_str());
    return request;
}

/************************************************************************/
/*                        DescribeCoverageRequest()                     */
/*                                                                      */
/************************************************************************/

CPLString WCSDataset201::DescribeCoverageRequest()
{
    CPLString request = CPLGetXMLValue(psService, "ServiceURL", "");
    request = CPLURLAddKVP(request, "SERVICE", "WCS");
    request = CPLURLAddKVP(request, "REQUEST", "DescribeCoverage");
    request = CPLURLAddKVP(request, "VERSION",
                           CPLGetXMLValue(psService, "Version", "2.0.1"));
    request = CPLURLAddKVP(request, "COVERAGEID",
                           CPLGetXMLValue(psService, "CoverageName", ""));
    request = CPLURLAddKVP(request, "FORMAT", "text/xml");
    CPLString extra = CPLGetXMLValue(psService, "Parameters", "");
    if (extra != "")
    {
        std::vector<CPLString> pairs = Split(extra, "&");
        for (unsigned int i = 0; i < pairs.size(); ++i)
        {
            std::vector<CPLString> pair = Split(pairs[i], "=");
            request = CPLURLAddKVP(request, pair[0], pair[1]);
        }
    }
    extra = CPLGetXMLValue(psService, "DescribeCoverageExtra", "");
    if (extra != "")
    {
        std::vector<CPLString> pairs = Split(extra, "&");
        for (unsigned int i = 0; i < pairs.size(); ++i)
        {
            std::vector<CPLString> pair = Split(pairs[i], "=");
            request = CPLURLAddKVP(request, pair[0], pair[1]);
        }
    }
    CPLDebug("WCS", "Requesting %s", request.c_str());
    return request;
}

/************************************************************************/
/*                             GridOffsets()                            */
/*                                                                      */
/************************************************************************/

bool WCSDataset201::GridOffsets(CPLXMLNode *grid, CPLString subtype,
                                bool swap_grid_axis,
                                std::vector<double> &origin,
                                std::vector<std::vector<double>> &offset,
                                std::vector<CPLString> axes, char ***metadata)
{
    // todo: use domain_index

    // origin position, center of cell
    CPLXMLNode *point = CPLGetXMLNode(grid, "origin.Point.pos");
    origin = Flist(
        Split(CPLGetXMLValue(point, nullptr, ""), " ", axis_order_swap), 0, 2);

    // offsets = coefficients of affine transformation from cell coords to
    // CRS coords, (1,2) and (4,5)

    if (subtype == "RectifiedGrid")
    {

        // for rectified grid the geo transform is from origin and offsetVectors
        int i = 0;
        for (CPLXMLNode *node = grid->psChild; node != nullptr;
             node = node->psNext)
        {
            if (node->eType != CXT_Element ||
                !EQUAL(node->pszValue, "offsetVector"))
            {
                continue;
            }
            offset.push_back(Flist(
                Split(CPLGetXMLValue(node, nullptr, ""), " ", axis_order_swap),
                0, 2));
            if (i == 1)
            {
                break;
            }
            i++;
        }
        if (offset.size() < 2)
        {
            // error or not?
            std::vector<double> x;
            x.push_back(1);
            x.push_back(0);
            std::vector<double> y;
            y.push_back(0);
            y.push_back(1);
            offset.push_back(x);
            offset.push_back(y);
        }
        // if axis_order_swap
        // the offset order should be swapped
        // Rasdaman does it
        // MapServer and GeoServer not
        if (swap_grid_axis)
        {
            std::vector<double> tmp = offset[0];
            offset[0] = offset[1];
            offset[1] = tmp;
        }
    }
    else
    {  // if (coverage_type == "ReferenceableGrid"(ByVector)) {

        // for vector referenceable grid the geo transform is from
        // offsetVector, coefficients, gridAxesSpanned, sequenceRule
        // in generalGridAxis.GeneralGridAxis
        for (CPLXMLNode *node = grid->psChild; node != nullptr;
             node = node->psNext)
        {
            CPLXMLNode *axis = CPLGetXMLNode(node, "GeneralGridAxis");
            if (!axis)
            {
                continue;
            }
            CPLString spanned = CPLGetXMLValue(axis, "gridAxesSpanned", "");
            int index = IndexOf(spanned, axes);
            if (index == -1)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "This is not a rectilinear grid(?).");
                return false;
            }
            CPLString coeffs = CPLGetXMLValue(axis, "coefficients", "");
            if (coeffs != "")
            {
                *metadata = CSLSetNameValue(
                    *metadata, CPLString().Printf("DIMENSION_%i_COEFFS", index),
                    coeffs);
            }
            CPLString order =
                CPLGetXMLValue(axis, "sequenceRule.axisOrder", "");
            CPLString rule = CPLGetXMLValue(axis, "sequenceRule", "");
            if (!(order == "+1" && rule == "Linear"))
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Grids with sequence rule '%s' and axis order '%s' "
                         "are not supported.",
                         rule.c_str(), order.c_str());
                return false;
            }
            CPLXMLNode *offset_node = CPLGetXMLNode(axis, "offsetVector");
            if (offset_node)
            {
                offset.push_back(
                    Flist(Split(CPLGetXMLValue(offset_node, nullptr, ""), " ",
                                axis_order_swap),
                          0, 2));
            }
            else
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Missing offset vector in grid axis.");
                return false;
            }
        }
        // todo: make sure offset order is the same as the axes order but see
        // above
    }
    if (origin.size() < 2 || offset.size() < 2)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Could not parse origin or offset vectors from grid.");
        return false;
    }
    return true;
}

/************************************************************************/
/*                             GetSubdataset()                          */
/*                                                                      */
/************************************************************************/

CPLString WCSDataset201::GetSubdataset(const CPLString &coverage)
{
    char **metadata = GDALPamDataset::GetMetadata("SUBDATASETS");
    CPLString subdataset;
    if (metadata != nullptr)
    {
        for (int i = 0; metadata[i] != nullptr; ++i)
        {
            char *key;
            CPLString url = CPLParseNameValue(metadata[i], &key);
            if (key != nullptr && strstr(key, "SUBDATASET_") &&
                strstr(key, "_NAME"))
            {
                if (coverage == CPLURLGetValue(url, "coverageId"))
                {
                    subdataset = key;
                    subdataset.erase(subdataset.find("_NAME"), 5);
                    CPLFree(key);
                    break;
                }
            }
            CPLFree(key);
        }
    }
    return subdataset;
}

/************************************************************************/
/*                             SetFormat()                              */
/*                                                                      */
/************************************************************************/

bool WCSDataset201::SetFormat(CPLXMLNode *coverage)
{
    // set the Format value in service,
    // unless it is set by the user
    CPLString format = CPLGetXMLValue(psService, "Format", "");

    // todo: check the value against list of supported formats?
    if (format != "")
    {
        return true;
    }

    /*      We will prefer anything that sounds like TIFF, otherwise        */
    /*      falling back to the first supported format.  Should we          */
    /*      consider preferring the nativeFormat if available?              */

    char **metadata = GDALPamDataset::GetMetadata(nullptr);
    const char *value =
        CSLFetchNameValue(metadata, "WCS_GLOBAL#formatSupported");
    if (value == nullptr)
    {
        format = CPLGetXMLValue(coverage, "ServiceParameters.nativeFormat", "");
    }
    else
    {
        std::vector<CPLString> format_list = Split(value, ",");
        for (unsigned j = 0; j < format_list.size(); ++j)
        {
            if (format_list[j].ifind("tiff") != std::string::npos)
            {
                format = format_list[j];
                break;
            }
        }
        if (format == "" && format_list.size() > 0)
        {
            format = format_list[0];
        }
    }
    if (format != "")
    {
        CPLSetXMLValue(psService, "Format", format);
        bServiceDirty = true;
        return true;
    }
    else
    {
        return false;
    }
}

/************************************************************************/
/*                         ParseGridFunction()                          */
/*                                                                      */
/************************************************************************/

bool WCSDataset201::ParseGridFunction(CPLXMLNode *coverage,
                                      std::vector<int> &axisOrder)
{
    CPLXMLNode *function =
        CPLGetXMLNode(coverage, "coverageFunction.GridFunction");
    if (function)
    {
        CPLString path = "sequenceRule";
        CPLString sequenceRule = CPLGetXMLValue(function, path, "");
        path += ".axisOrder";
        axisOrder = Ilist(Split(CPLGetXMLValue(function, path, ""), " "));
        // for now require simple
        if (sequenceRule != "Linear")
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Can't handle '%s' coverages.", sequenceRule.c_str());
            return false;
        }
    }
    return true;
}

/************************************************************************/
/*                             ParseRange()                             */
/*                                                                      */
/************************************************************************/

int WCSDataset201::ParseRange(CPLXMLNode *coverage,
                              const CPLString &range_subset, char ***metadata)
{
    int fields = 0;
    // Default is to include all (types permitting?)
    // Can also be controlled with Range parameter

    // The contents of a rangeType is a swe:DataRecord
    CPLString path = "rangeType.DataRecord";
    CPLXMLNode *record = CPLGetXMLNode(coverage, path);
    if (!record)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Attributes are not defined in a DataRecord, giving up.");
        return 0;
    }

    // mapserver does not like field names, it wants indexes
    // so we should be able to give those

    // if Range is set remove those not in it
    std::vector<CPLString> range = Split(range_subset, ",");
    // todo: add check for range subsetting profile existence in server metadata
    // here
    unsigned int range_index = 0;  // index for reading from range
    bool in_band_range = false;

    unsigned int field_index = 1;
    CPLString field_name;
    std::vector<CPLString> nodata_array;

    for (CPLXMLNode *field = record->psChild; field != nullptr;
         field = field->psNext)
    {
        if (field->eType != CXT_Element || !EQUAL(field->pszValue, "field"))
        {
            continue;
        }
        CPLString fname = CPLGetXMLValue(field, "name", "");
        bool include = true;

        if (range.size() > 0)
        {
            include = false;
            if (range_index < range.size())
            {
                CPLString current_range = range[range_index];
                CPLString fname_test;

                if (atoi(current_range) != 0)
                {
                    fname_test = CPLString().Printf("%i", field_index);
                }
                else
                {
                    fname_test = fname;
                }

                if (current_range == "*")
                {
                    include = true;
                }
                else if (current_range == fname_test)
                {
                    include = true;
                    range_index += 1;
                }
                else if (current_range.find(fname_test + ":") !=
                         std::string::npos)
                {
                    include = true;
                    in_band_range = true;
                }
                else if (current_range.find(":" + fname_test) !=
                         std::string::npos)
                {
                    include = true;
                    in_band_range = false;
                    range_index += 1;
                }
                else if (in_band_range)
                {
                    include = true;
                }
            }
        }

        if (include)
        {
            CPLString key;
            key.Printf("FIELD_%i_", field_index);
            *metadata =
                CSLSetNameValue(*metadata, (key + "NAME").c_str(), fname);

            CPLString nodata =
                CPLGetXMLValue(field, "Quantity.nilValues.NilValue", "");
            if (nodata != "")
            {
                *metadata = CSLSetNameValue(*metadata, (key + "NODATA").c_str(),
                                            nodata);
            }

            CPLString descr = CPLGetXMLValue(field, "Quantity.description", "");
            if (descr != "")
            {
                *metadata =
                    CSLSetNameValue(*metadata, (key + "DESCR").c_str(), descr);
            }

            path = "Quantity.constraint.AllowedValues.interval";
            CPLString interval = CPLGetXMLValue(field, path, "");
            if (interval != "")
            {
                *metadata = CSLSetNameValue(
                    *metadata, (key + "INTERVAL").c_str(), interval);
            }

            if (field_name == "")
            {
                field_name = fname;
            }
            nodata_array.push_back(nodata);
            fields += 1;
        }

        field_index += 1;
    }

    if (fields == 0)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "No data fields found (bad Range?).");
    }
    else
    {
        // todo: default to the first one?
        bServiceDirty =
            CPLUpdateXML(psService, "NoDataValue", Join(nodata_array, ",")) ||
            bServiceDirty;
    }

    return fields;
}

/************************************************************************/
/*                          ExtractGridInfo()                           */
/*                                                                      */
/*      Collect info about grid from describe coverage for WCS 2.0.     */
/*                                                                      */
/************************************************************************/

bool WCSDataset201::ExtractGridInfo()
{
    // this is for checking what's in service and for filling in empty slots in
    // it if the service file can be considered ready for use, this could be
    // skipped

    CPLXMLNode *coverage = CPLGetXMLNode(psService, "CoverageDescription");

    if (coverage == nullptr)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "CoverageDescription missing from service.");
        return false;
    }

    CPLString subtype = CoverageSubtype(coverage);

    // get CRS from boundedBy.Envelope and set the native flag to true
    // below we may set the CRS again but that won't be native (however, non
    // native CRS is not yet supported) also axis order swap is set
    CPLString path = "boundedBy.Envelope";
    CPLXMLNode *envelope = CPLGetXMLNode(coverage, path);
    if (envelope == nullptr)
    {
        path = "boundedBy.EnvelopeWithTimePeriod";
        envelope = CPLGetXMLNode(coverage, path);
        if (envelope == nullptr)
        {
            CPLError(CE_Failure, CPLE_AppDefined, "Missing boundedBy.Envelope");
            return false;
        }
    }
    std::vector<CPLString> bbox = ParseBoundingBox(envelope);
    if (!SetCRS(ParseCRS(envelope), true) || bbox.size() < 2)
    {
        return false;
    }

    // has the user set the domain?
    std::vector<CPLString> domain =
        Split(CPLGetXMLValue(psService, "Domain", ""), ",");

    // names of axes
    std::vector<CPLString> axes =
        Split(CPLGetXMLValue(coverage, (path + ".axisLabels").c_str(), ""), " ",
              axis_order_swap);
    std::vector<CPLString> uoms =
        Split(CPLGetXMLValue(coverage, (path + ".uomLabels").c_str(), ""), " ",
              axis_order_swap);

    if (axes.size() < 2)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "The coverage has less than 2 dimensions or no axisLabels.");
        return false;
    }

    std::vector<int> domain_indexes = IndexOf(domain, axes);
    if (Contains(domain_indexes, -1))
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Axis in given domain does not exist in coverage.");
        return false;
    }
    if (domain_indexes.size() == 0)
    {  // default is the first two
        domain_indexes.push_back(0);
        domain_indexes.push_back(1);
    }
    if (domain.size() == 0)
    {
        domain.push_back(axes[0]);
        domain.push_back(axes[1]);
        CPLSetXMLValue(psService, "Domain", Join(domain, ","));
        bServiceDirty = true;
    }

    // GridFunction (is optional)
    // We support only linear grid functions.
    // axisOrder determines how data is arranged in the grid <order><axis
    // number> specifically: +2 +1 => swap grid envelope and the order of the
    // offsets
    std::vector<int> axisOrder;
    if (!ParseGridFunction(coverage, axisOrder))
    {
        return false;
    }

    const char *md_domain = "";
    char **metadata = CSLDuplicate(
        GetMetadata(md_domain));  // coverage metadata to be added/updated

    metadata = CSLSetNameValue(metadata, "DOMAIN", Join(domain, ","));

    // add coverage metadata: GeoServer TimeDomain

    CPLXMLNode *timedomain =
        CPLGetXMLNode(coverage, "metadata.Extension.TimeDomain");
    if (timedomain)
    {
        std::vector<CPLString> timePositions;
        // "//timePosition"
        for (CPLXMLNode *node = timedomain->psChild; node != nullptr;
             node = node->psNext)
        {
            if (node->eType != CXT_Element ||
                strcmp(node->pszValue, "TimeInstant") != 0)
            {
                continue;
            }
            for (CPLXMLNode *node2 = node->psChild; node2 != nullptr;
                 node2 = node2->psNext)
            {
                if (node2->eType != CXT_Element ||
                    strcmp(node2->pszValue, "timePosition") != 0)
                {
                    continue;
                }
                timePositions.push_back(CPLGetXMLValue(node2, "", ""));
            }
        }
        metadata =
            CSLSetNameValue(metadata, "TimeDomain", Join(timePositions, ","));
    }

    // dimension metadata

    std::vector<CPLString> slow = Split(bbox[0], " ", axis_order_swap);
    std::vector<CPLString> shigh = Split(bbox[1], " ", axis_order_swap);
    bServiceDirty =
        CPLUpdateXML(psService, "Low", Join(slow, ",")) || bServiceDirty;
    bServiceDirty =
        CPLUpdateXML(psService, "High", Join(shigh, ",")) || bServiceDirty;
    if (slow.size() < 2 || shigh.size() < 2)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "The coverage has less than 2 dimensions.");
        CSLDestroy(metadata);
        return false;
    }
    // todo: if our x,y domain is not the first two? use domain_indexes?
    std::vector<double> low = Flist(slow, 0, 2);
    std::vector<double> high = Flist(shigh, 0, 2);
    std::vector<double> env;
    env.insert(env.end(), low.begin(), low.begin() + 2);
    env.insert(env.end(), high.begin(), high.begin() + 2);

    for (unsigned int i = 0; i < axes.size(); ++i)
    {
        CPLString key;
        key.Printf("DIMENSION_%i_", i);
        metadata = CSLSetNameValue(metadata, (key + "AXIS").c_str(), axes[i]);
        if (i < uoms.size())
        {
            metadata =
                CSLSetNameValue(metadata, (key + "UOM").c_str(), uoms[i]);
        }
        if (i < 2)
        {
            metadata = CSLSetNameValue(
                metadata, (key + "INTERVAL").c_str(),
                CPLString().Printf("%.15g,%.15g", low[i], high[i]));
        }
        else if (i < slow.size() && i < shigh.size())
        {
            metadata = CSLSetNameValue(
                metadata, (key + "INTERVAL").c_str(),
                CPLString().Printf("%s,%s", slow[i].c_str(), shigh[i].c_str()));
        }
        else if (i < bbox.size())
        {
            metadata = CSLSetNameValue(metadata, (key + "INTERVAL").c_str(),
                                       bbox[i].c_str());
        }
    }

    // domainSet
    // requirement 23: the srsName here _shall_ be the same as in boundedBy
    // => we ignore it
    // the CRS of this dataset is from boundedBy (unless it is overridden)
    // this is the size of this dataset
    // this gives the geotransform of this dataset (unless there is CRS
    // override)

    CPLXMLNode *grid = GetGridNode(coverage, subtype);
    if (!grid)
    {
        CSLDestroy(metadata);
        return false;
    }

    //
    bool swap_grid_axis = false;
    if (axisOrder.size() >= 2 && axisOrder[domain_indexes[0]] == 2 &&
        axisOrder[domain_indexes[1]] == 1)
    {
        swap_grid_axis = !CPLGetXMLBoolean(psService, "NoGridAxisSwap");
    }
    path = "limits.GridEnvelope";
    std::vector<std::vector<int>> size =
        ParseGridEnvelope(CPLGetXMLNode(grid, path), swap_grid_axis);
    std::vector<int> grid_size;
    if (size.size() < 2)
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Can't parse the grid envelope.");
        CSLDestroy(metadata);
        return false;
    }

    grid_size.push_back(size[1][domain_indexes[0]] -
                        size[0][domain_indexes[0]] + 1);
    grid_size.push_back(size[1][domain_indexes[1]] -
                        size[0][domain_indexes[1]] + 1);

    path = "axisLabels";
    bool swap_grid_axis_labels =
        swap_grid_axis || CPLGetXMLBoolean(psService, "GridAxisLabelSwap");
    std::vector<CPLString> grid_axes =
        Split(CPLGetXMLValue(grid, path, ""), " ", swap_grid_axis_labels);
    // autocorrect MapServer thing
    if (grid_axes.size() >= 2 && grid_axes[0] == "lat" &&
        grid_axes[1] == "long")
    {
        grid_axes[0] = "long";
        grid_axes[1] = "lat";
    }
    bServiceDirty = CPLUpdateXML(psService, "GridAxes", Join(grid_axes, ",")) ||
                    bServiceDirty;

    std::vector<double> origin;
    std::vector<std::vector<double>> offsets;
    if (!GridOffsets(grid, subtype, swap_grid_axis, origin, offsets, axes,
                     &metadata))
    {
        CSLDestroy(metadata);
        return false;
    }

    SetGeometry(grid_size, origin, offsets);

    // subsetting and dimension to bands
    std::vector<CPLString> dimensions;
    CPLString range;
    std::vector<std::vector<CPLString>> others;
    ParseParameters(psService, dimensions, range, others);

    // it is ok to have trimming or even slicing for x/y, it just affects our
    // bounding box but that is a todo item todo: BoundGeometry(domain_trim) if
    // domain_trim.size() > 0
    std::vector<std::vector<double>> domain_trim;

    // are all dimensions that are not x/y domain sliced?
    // if not, bands can't be defined, see below
    bool dimensions_are_ok = true;
    for (unsigned int i = 0; i < axes.size(); ++i)
    {
        std::vector<CPLString> params;
        for (unsigned int j = 0; j < dimensions.size(); ++j)
        {
            if (dimensions[j].find(axes[i] + "(") != std::string::npos)
            {
                params = Split(FromParenthesis(dimensions[j]), ",");
                break;
            }
        }
        int domain_index = IndexOf(axes[i], domain);
        if (domain_index != -1)
        {
            std::vector<double> trim = Flist(params, 0, 2);
            domain_trim.push_back(trim);
            continue;
        }
        // size == 1 => sliced
        if (params.size() != 1)
        {
            dimensions_are_ok = false;
        }
    }
    // todo: add metadata: note: no bands, you need to subset to get data

    // check for CRS override
    CPLString crs = CPLGetXMLValue(psService, "SRS", "");
    if (crs != "" && crs != osCRS)
    {
        if (!SetCRS(crs, false))
        {
            CSLDestroy(metadata);
            return false;
        }
        // todo: support CRS override, it requires warping the grid to the new
        // CRS SetGeometry(grid_size, origin, offsets);
        CPLError(CE_Failure, CPLE_AppDefined,
                 "CRS override not yet supported.");
        CSLDestroy(metadata);
        return false;
    }

    // todo: ElevationDomain, DimensionDomain

    // rangeType

    // get the field metadata
    // get the count of fields
    // if Range is set in service that may limit the fields
    int fields = ParseRange(coverage, range, &metadata);
    // if fields is 0 an error message has been emitted
    // but we let this go on since the user may be experimenting
    // and she wants to see the resulting metadata and not just an error message
    // situation is ~the same when bands == 0 when we exit here

    // todo: do this only if metadata is dirty
    this->SetMetadata(metadata, md_domain);
    CSLDestroy(metadata);
    TrySaveXML();

    // determine the band count
    int bands = 0;
    if (dimensions_are_ok)
    {
        bands = fields;
    }
    bServiceDirty =
        CPLUpdateXML(psService, "BandCount", CPLString().Printf("%d", bands)) ||
        bServiceDirty;

    // set the Format value in service, unless it is set
    // by the user, either through direct edit or options
    if (!SetFormat(coverage))
    {
        // all attempts to find a format have failed...
        CPLError(CE_Failure, CPLE_AppDefined,
                 "All attempts to find a format have failed, giving up.");
        return false;
    }

    return true;
}