Fix: WCS DescribeCoverage does not show dimensions for a multidim layer

CMakeLists.txt

@@ -50,6 +50,7 @@ add_executable(geojsondump ./adagucserverEC/geojsondump.cpp)
 
 add_test(testhclasses hclasses/testhclasses)
 add_test(testadagucserver adagucserverEC/testadagucserver)
+add_test(testtimeutils adagucserverEC/testtimeutils)
 
 target_link_libraries(adagucserver adagucserverEC hclasses CCDFDataModel)
 target_link_libraries(h5ncdump adagucserverEC hclasses CCDFDataModel)

adagucserverEC/CMakeLists.txt

@@ -196,6 +196,9 @@ add_library(
     Types/ProjectionStore.h
     Types/ProjectionStore.cpp 
     testadagucserver.cpp
+    testtimeutils.cpp
+    timeutils.cpp
+    timeutils.h
 )
 
 target_include_directories(adagucserverEC PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${Cairo_INCLUDE_DIRS} ${FREETYPE_INCLUDE_DIRS} ${PostgreSQL_INCLUDE_DIRS} ${GDAL_INCLUDE_DIRS} ${PROJ_INCLUDE_DIR})
@@ -204,3 +207,7 @@ target_link_libraries(adagucserverEC CCDFDataModel hclasses ${NetCDF_LIBRARIES}
 # Build unit test executable
 add_executable(testadagucserver testadagucserver.cpp)
 target_link_libraries(testadagucserver PRIVATE adagucserverEC CCDFDataModel hclasses CppUnitLite)
+
+# Build unit tests for timeutils
+add_executable(testtimeutils testtimeutils.cpp)
+target_link_libraries(testtimeutils PRIVATE adagucserverEC hclasses CppUnitLite)

adagucserverEC/CXMLGen.cpp

@@ -25,10 +25,12 @@
 
 #include <algorithm>
 #include <vector>
+#include <sstream>
 #include <string>
 #include "CXMLGen.h"
 #include "CDBFactory.h"
 #include "LayerTypeLiveUpdate/LayerTypeLiveUpdate.h"
+#include "timeutils.h"
 #include <json_adaguc.h>
 #include "utils/LayerMetadataStore.h"
 #include "utils/XMLGenUtils.h"
@@ -41,6 +43,35 @@
 const char *CXMLGen::className = "CXMLGen";
 int CXMLGen::WCSDescribeCoverage(CServerParams *srvParam, CT::string *XMLDocument) { return OGCGetCapabilities(srvParam, XMLDocument); }
 
+// Function to parse a string to double if numeric
+double parseNumeric(std::string const &str, bool &isNumeric) {
+  auto result = double();
+  auto i = std::istringstream(str);
+  i >> result;
+  isNumeric = !i.fail() && i.eof();
+  return result;
+}
+
+// Sort values that can either be numeric of a string
+bool multiTypeSort(const CT::string &a, const CT::string &b) {
+  // Try to convert strings to numbers
+  float aNum, bNum;
+  bool isANum, isBNum;
+
+  isANum = false;
+  aNum = parseNumeric(a.c_str(), isANum);
+
+  isBNum = false;
+  bNum = parseNumeric(b.c_str(), isBNum);
+
+  // Do numerical comparison or alphabetical comparison according to type
+  if (isANum && isBNum) {
+    return aNum < bNum;
+  } else {
+    return a < b;
+  }
+}
+
 const MetadataLayer *getFirstLayerWithoutError(std::vector<MetadataLayer *> *metadataLayerList) {
   if (metadataLayerList->size() == 0) {
     return nullptr;
@@ -775,6 +806,79 @@ int CXMLGen::getWCS_1_0_0_Capabilities(CT::string *XMLDoc, std::vector<MetadataL
   return 0;
 }
 
+void generateRangeSet(CT::string *XMLDoc, MetadataLayer *layer) {
+  /*
+  From the documentation:
+  The optional and repeatable axisDescription/AxisDescription element is for compound observations.
+  It describes an additional parameter (that is, an independent variable besides space and time),
+  and the valid values of this parameter, which GetCoverage requests can use to select subsets of a
+  coverage offering.
+  */
+  if (!layer->layerMetadata.dimList.size()) {
+    return;
+  }
+  XMLDoc->concat("    <rangeSet>\n"
+                 "      <RangeSet>\n"
+                 "        <name>dimensions</name>\n"
+                 "        <label>dimensions</label>\n");
+  // Dims
+  for (size_t d = 0; d < layer->layerMetadata.dimList.size(); d++) {
+    LayerMetadataDim *dim = &layer->layerMetadata.dimList[d];
+    CT::string min, max, duration;
+    CT::string *valueSplit;
+    std::vector<CT::string> valuesVector;
+
+    // Case of min/max(/duration), for time dimension
+    valueSplit = dim->values.splitToArray("/");
+    if (valueSplit->count >= 2) {
+      min = valueSplit[0];
+      max = valueSplit[1];
+      // Third value is the interval duration
+      if (valueSplit->count == 3) {
+        duration = valueSplit[2];
+      }
+    } else {
+      // General case of a list of values (of any type)
+      valueSplit = dim->values.splitToArray(",");
+      valuesVector = std::vector<CT::string>(valueSplit, valueSplit + valueSplit->count);
+      std::sort(valuesVector.begin(), valuesVector.end(), multiTypeSort);
+      min = valuesVector[0];
+      max = valuesVector.back();
+    }
+
+    XMLDoc->printconcat("        <axisDescription>\n"
+                        "          <AxisDescription>\n"
+                        "            <name>%s</name>\n"
+                        "            <label>%s</label>\n",
+                        dim->cdfName.c_str(), dim->cdfName.c_str());
+    if (valueSplit->count >= 2) {
+      XMLDoc->printconcat("            <values>\n"
+                          "              <interval>\n"
+                          "                <min>%s</min>\n"
+                          "                <max>%s</max>\n",
+                          min.c_str(), max.c_str());
+      // Precalculate the interval in the case of time (no interval if fewer than 4 values)
+      if ((dim->cdfName.indexOf("time") != -1) && duration.length() > 0) {
+        XMLDoc->printconcat("                <res>%s</res>\n", duration.c_str()); // .c_str());
+      }
+      XMLDoc->printconcat("              </interval>\n");
+      // Print all possible values if there is a relatively small number, for other dimensions
+      if ((valueSplit->count <= 100) && (dim->cdfName.indexOf("time") == -1)) {
+        for (size_t i = 0; i < valueSplit->count; i++) {
+          XMLDoc->printconcat("              <singleValue>%s</singleValue>\n", valuesVector[i].c_str());
+        }
+      }
+      XMLDoc->printconcat("            </values>\n");
+    }
+    XMLDoc->printconcat("          </AxisDescription>\n"
+                        "        </axisDescription>\n");
+    delete[] valueSplit;
+  }
+
+  XMLDoc->concat("      </RangeSet>\n"
+                 "    </rangeSet>\n");
+}
+
 int CXMLGen::getWCS_1_0_0_DescribeCoverage(CT::string *XMLDoc, std::vector<MetadataLayer *> *metadataLayerList) {
 
   XMLDoc->copy("<?xml version='1.0' encoding=\"ISO-8859-1\" ?>\n"
@@ -902,22 +1006,9 @@ int CXMLGen::getWCS_1_0_0_DescribeCoverage(CT::string *XMLDoc, std::vector<Metad
                 }
                 XMLDoc->concat("      </temporalDomain>\n");
               }
-              XMLDoc->concat("    </domainSet>\n"
-                             "    <rangeSet>\n"
-                             "      <RangeSet>\n"
-                             "        <name>bands</name>\n"
-                             "        <label>bands</label>\n"
-                             "        <axisDescription>\n"
-                             "          <AxisDescription>\n"
-                             "            <name>bands</name>\n"
-                             "            <label>Bands/Channels/Samples</label>\n"
-                             "            <values>\n"
-                             "              <singleValue>1</singleValue>\n"
-                             "            </values>\n"
-                             "          </AxisDescription>\n"
-                             "        </axisDescription>\n"
-                             "      </RangeSet>\n"
-                             "    </rangeSet>\n");
+              XMLDoc->concat("    </domainSet>\n");
+              // Generate the XML code for RangeSet, including dimensions (AxisDescriptions)
+              generateRangeSet(XMLDoc, layer);
               // Supported CRSs
               XMLDoc->concat("    <supportedCRSs>\n");
 

adagucserverEC/testtimeutils.cpp

@@ -0,0 +1,103 @@
+/******************************************************************************
+ *
+ * Copyright 2024, Royal Netherlands Meteorological Institute (KNMI)
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ ******************************************************************************/
+
+#include "CDebugger.h"
+#include "CppUnitLite/TestHarness.h"
+#include "CTString.h"
+#include "timeutils.h"
+
+TEST(CalculateTimeInterval, TimeUtils) {
+  CTime::Date start = {};
+  start.minute = 30;
+  start.hour = 15;
+  start.day = 29;
+  start.month = 8;
+  start.year = 2024;
+  start.second = 0.0;
+
+  CTime::Date end = {};
+  end.second = 45.0;
+  end.minute = 45;
+  end.hour = 16;
+  end.day = 29;
+  end.month = 8;
+  end.year = 2025;
+
+  TimeInterval expected = {1, 0, 0, 1, 15, 45};
+  TimeInterval result = calculateTimeInterval(start, end);
+
+  CHECK(result.years == expected.years);
+  CHECK(result.months == expected.months);
+  CHECK(result.days == expected.days);
+  CHECK(result.hours == expected.hours);
+  CHECK(result.minutes == expected.minutes);
+  CHECK(result.seconds == expected.seconds);
+}
+
+TEST(ToISO8601Interval, TimeUtils) {
+  TimeInterval interval = {1, 2, 10, 5, 30, 45};
+  CT::string expected("P1Y2M10DT5H30M45S");
+
+  CT::string result = toISO8601Interval(interval);
+  CHECK(result == expected);
+}
+
+TEST(ToSeconds, TimeUtils) {
+  TimeInterval interval = {0, 1, 0, 1, 0, 20}; // 1 month, 1 hour, 20 seconds
+  long expected = 2595620;
+
+  long result = interval.toSeconds();
+  CHECK(result == expected);
+}
+
+TEST(EstimateISO8601Duration, TimeUtils) {
+  // Most straightforward case
+  std::vector<CT::string> timestamps = {"2024-07-29T15:30:45Z", "2024-07-29T16:30:45Z", "2024-07-29T17:30:45Z", "2024-07-29T18:30:45Z"};
+  CT::string expected("PT1H");
+  CT::string result = estimateISO8601Duration(timestamps);
+  std::cout << "Result: " << result << "\n";
+  CHECK(result == expected);
+
+  // Case with one timestep missing
+  std::vector<CT::string> timestamps_with_gap = {"2024-07-29T15:30:45Z", "2024-07-29T16:30:45Z", "2024-07-29T18:30:45Z", "2024-07-29T19:30:45Z",
+                                                 "2024-07-29T20:30:45Z", "2024-07-29T21:30:45Z", "2024-07-29T22:30:45Z", "2024-07-29T23:30:45Z"};
+  CT::string expected_with_gap("PT1H");
+  CT::string result_with_gap = estimateISO8601Duration(timestamps_with_gap);
+  CHECK(result_with_gap == expected_with_gap);
+
+  // Case in February of a leap year
+  std::vector<CT::string> leap_year_timestamps = {"2024-02-28T23:30:45Z", "2024-02-29T00:30:45Z", "2024-02-29T01:30:45Z", "2024-02-29T02:30:45Z"};
+  CT::string expected_leap_year("PT1H");
+  CT::string result_leap_year = estimateISO8601Duration(leap_year_timestamps);
+  CHECK(result_leap_year == expected_leap_year);
+
+  // Case of a complex interval (5 days and 10 hours)
+  std::vector<CT::string> complex_interval_timestamps = {"2024-07-01T08:00:00Z", "2024-07-06T18:00:00Z", "2024-07-12T04:00:00Z", "2024-07-17T14:00:00Z", "2024-07-23T00:00:00Z"};
+  CT::string expected_complex_interval("P5DT10H");
+  CT::string result_complex_interval = estimateISO8601Duration(complex_interval_timestamps);
+  CHECK(result_complex_interval == expected_complex_interval);
+}
+
+int main() {
+  TestResult tr;
+  TestRegistry::runAllTests(tr);
+  if (tr.failureCount != 0) {
+    return 1;
+  }
+  return 0;
+}