REFAC: make otb::SarSensorModel independant of otb::ImageMetadata and use MJD...

REFAC: make otb::SarSensorModel independant of otb::ImageMetadata and use MJD stored as double (adapted from ossim plugins) instead of boost datetime (which has insufficient microsecond precision by default)

Modules/Adapters/BoostAdapters/CMakeLists.txt

@@ -19,7 +19,4 @@
 #
 
 project(OTBBoostAdapters)
-
-set(OTBBoostAdapters_LIBRARIES OTBBoostAdapters)
-
 otb_module_impl()

Modules/Adapters/BoostAdapters/otb-module.cmake

@@ -21,7 +21,6 @@
 set(DOCUMENTATION "Adapters for the Boost library.")
 
 otb_module(OTBBoostAdapters 
-ENABLE_SHARED
   DEPENDS
     OTBBoost
   DESCRIPTION

Modules/Adapters/BoostAdapters/src/CMakeLists.txt

-#
-# Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
-#
-# This file is part of Orfeo Toolbox
-#
-#     https://www.orfeo-toolbox.org/
-#
-# 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.
-#
-
-set(OTBBoostAdapters_SRC
-  otbDateTime.cxx
-  )
-
-add_library(OTBBoostAdapters ${OTBBoostAdapters_SRC})
-target_link_libraries(OTBBoostAdapters 
-  ${OTBBoost_LIBRARIES}
-  )
-
-otb_module_target(OTBBoostAdapters)

Modules/Adapters/BoostAdapters/src/otbDateTime.cxx

-/*
- * Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
- *
- * This file is part of Orfeo Toolbox
- *
- *     https://www.orfeo-toolbox.org/
- *
- * 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 "otbDateTime.h"
-
-namespace boost 
-{ 
-namespace posix_time 
-{
-
-time_duration abs(time_duration d)
-{
-  if(d.is_negative())
-    d = d.invert_sign();
-  return d;
-}
-
-}
-}
\ No newline at end of file

Modules/Adapters/BoostAdapters/include/otbDateTime.h → Modules/Core/Metadata/include/otbDateTime.h

@@ -21,6 +21,10 @@
 #ifndef ossimDateTime_h
 #define ossimDateTime_h
 
+#define OTB_USE_BOOST_TIME 0
+
+#if OTB_USE_BOOST_TIME
+
 #include <boost/date_time/posix_time/posix_time.hpp>
 
 #include <boost/config.hpp>
@@ -32,10 +36,10 @@
 // Hence this new class injected into boost namespace to emulate a duration
 // with precision behind the microsecond.
 namespace boost { namespace posix_time {
-   class precise_duration;
-   double ratio_(precise_duration const& lhs, precise_duration const& rhs);
+  class precise_duration;
+  double ratio_(precise_duration const& lhs, precise_duration const& rhs);
 
-   class precise_duration
+  class precise_duration
       : private boostAdapter::addable<precise_duration>
       , private boostAdapter::substractable<precise_duration>
       , private boostAdapter::streamable<precise_duration>
@@ -139,9 +143,194 @@ namespace otb
 {
 namespace MetaData
 {
-  using TimeType = boost::posix_time::ptime;
-  using DurationType = boost::posix_time::precise_duration;
+using TimeType = boost::posix_time::ptime;
+using DurationType = boost::posix_time::precise_duration;
+
+DurationType seconds(double);
+
+TimeType ReadFormattedDate(const std::string & dateStr, const std::string & format = "%Y-%m-%dT%H:%M:%S%F");
+
+}
+}
+
+#else //OTB_USE_BOOST_TIME
+
+#include "otbOperatorUtilities.h"
+#include <cassert>
+#include <iomanip>
+
+namespace otb
+{
+namespace MetaData
+{
+namespace details
+{
+class DayFrac
+{
+  public:
+    typedef double scalar_type;
+    // DayFrac(DayFrac const&) = default;
+    // DayFrac(DayFrac &&) = default;
+         // DayFrac& operator=(DayFrac const&) = default;
+         // DayFrac& operator=(DayFrac &&) = default;
+    double as_day_frac() const { return m_day_frac; }
+
+    std::ostream & display(std::ostream & os) const { return os << m_day_frac; }
+    std::istream & read   (std::istream & is)       { return is >> m_day_frac; }
+
+      protected:
+         /**@name Construction/destruction
+         */
+         //@{
+         /** Initialization constructor.
+         */
+         explicit DayFrac() = default;
+         explicit DayFrac(double day_frac) : m_day_frac(day_frac) {}
+         /** Protected destructor.
+         */
+         ~DayFrac() = default;
+         //@}
+
+         /**@name Operations
+         */
+         //@{
+         void add(DayFrac const& rhs) { m_day_frac += rhs.m_day_frac; }
+         void sub(DayFrac const& rhs) { m_day_frac -= rhs.m_day_frac; }
+         void mult(scalar_type coeff) { m_day_frac *= coeff; }
+         void div(scalar_type coeff)  { assert(coeff && "Cannot divide by 0"); m_day_frac /= coeff; }
+         template <typename V> friend scalar_type ratio_(V const& lhs, V const& rhs)
+         { return lhs.as_day_frac() / rhs.as_day_frac(); }
+
+         template <typename U, typename V> friend U& operator+=(U & u, V const& v) {
+            u.add(v);
+            return u;
+         }
+         template <typename U, typename V> friend U& operator-=(U & u, V const& v) {
+            u.sub(v);
+            return u;
+         }
+
+         template <typename U, typename V> static U diff(V const& lhs, V const& rhs) {
+            U const res(lhs.as_day_frac() - rhs.as_day_frac());
+            return res;
+         }
+
+         template <typename U> friend U& operator*=(U & u, scalar_type const& v) {
+            u.mult(v);
+            return u;
+         }
+         template <typename U> friend U& operator/=(U & u, scalar_type const& v) {
+            u.div(v);
+            return u;
+         }
+
+         template <typename T> friend bool operator<(T const& lhs, T const& rhs) {
+            return lhs.as_day_frac() < rhs.as_day_frac();
+         }
+         template <typename T> friend bool operator==(T const& lhs, T const& rhs) {
+            return lhs.as_day_frac() == rhs.as_day_frac();
+         }
+         //@}
+      private:
+         double m_day_frac;
+      };
+   }
+
+   /**
+    * Duration abstraction.
+    *
+    * Values of this class represent time interval.
+    *
+    * <p><b>Semantics</b><br>
+    * <li> Value, mathematical: it's relative position
+    * <li> Time interval
+    *
+    * @see \c std::duration<>
+    */
+   class Duration
+      : public details::DayFrac
+      , private boostAdapter::addable<Duration>
+      , private boostAdapter::substractable<Duration>
+      , private boostAdapter::streamable<Duration>
+      , private boostAdapter::multipliable2<Duration, double>
+      , private boostAdapter::dividable<Duration, details::DayFrac::scalar_type>
+      , private boostAdapter::equality_comparable<Duration>
+      , private boostAdapter::less_than_comparable<Duration>
+      {
+      public:
+         typedef details::DayFrac::scalar_type scalar_type;
+
+         /**@name Construction/destruction
+         */
+         //@{
+         /** Initialization constructor.
+         */
+         Duration() = default;
+         explicit Duration(double day_frac)
+            : details::DayFrac(day_frac) {}
+         //@}
+
+         double total_seconds() const {
+            return as_day_frac() * 24 * 60 * 60;
+         }
+         double total_microseconds() const {
+            return total_seconds() * 1000 * 1000;
+         }
+         bool is_negative() const { return as_day_frac() < 0.0; }
+         Duration invert_sign() { return Duration(- as_day_frac()); }
+         friend Duration abs(Duration const& d) { return Duration(std::abs(d.as_day_frac())); }
+      };
+
+   /**
+    * Modified JulianDate abstraction.
+    *
+    * Objects of this class represent points in time.
+    *
+    * <p><b>Semantics</b><br>
+    * <li> Value, mathematical: it's an absolute position
+    * <li> Point in time
+    * @see \c std::time_point<>
+    */
+   class ModifiedJulianDate
+      : public details::DayFrac
+      , private boostAdapter::addable<ModifiedJulianDate, Duration>
+      , private boostAdapter::substractable<ModifiedJulianDate, Duration>
+      , private boostAdapter::substractable_asym<Duration, ModifiedJulianDate>
+      , private boostAdapter::streamable<ModifiedJulianDate>
+      , private boostAdapter::equality_comparable<ModifiedJulianDate>
+      , private boostAdapter::less_than_comparable<ModifiedJulianDate>
+      {
+      public:
+         typedef details::DayFrac::scalar_type scalar_type;
+
+         /**@name Construction/destruction
+         */
+         //@{
+         /** Initialization constructor.
+         */
+         ModifiedJulianDate() = default;
+         explicit ModifiedJulianDate(double day_frac)
+            : details::DayFrac(day_frac) {}
+         //@}
+         using details::DayFrac::diff;
+      };
+
+   ModifiedJulianDate toModifiedJulianDate(std::string const& utcTimeString);
+   inline Duration microseconds(double us) {
+      return Duration(us / (24ULL * 60 * 60 * 1000 * 1000));
+  }
+  inline Duration seconds(double us) {
+    return Duration(us / (24ULL * 60 * 60));
+  }
+
+  ModifiedJulianDate ReadFormattedDate(const std::string & dateStr, const std::string & format = "%Y-%m-%dT%H:%M:%S%F");
+
+  using TimeType = ModifiedJulianDate;
+  using DurationType = Duration;
 }
 }
 
 #endif
+
+
+#endif

Modules/Core/Metadata/include/otbSentinel1ImageMetadataInterface.h

@@ -100,6 +100,12 @@ public:
   bool HasCalibrationLookupDataFlag(const MetadataSupplierInterface&) const override;
   bool CreateCalibrationLookupData(SARCalib&, const ImageMetadata&, const MetadataSupplierInterface&, const bool) const override;
 
+  /* Read the Sar parameters and GCPs from the annotation xml file, this method can be used 
+    to instantiate a otb::SarModel */
+  void ReadSarParamAndGCPs(const XMLMetadataSupplier &,
+                           SARParam &,
+                           Projection::GCPParam &);
+
   void ParseGdal(ImageMetadata &) override;
 
   void ParseGeom(ImageMetadata &) override;
@@ -133,10 +139,6 @@ protected:
                                                                         const std::string & rg0_path,
                                                                         const std::string & coeffs_path) const;
 
-  /* Fetch azimuth and range times corresponding to GCPs */
-  std::unordered_map<std::string, GCPTime> GetGCPTimes(const XMLMetadataSupplier & xmlMS,
-                                                       const Projection::GCPParam & gcps) const;
-
   /* Compute the mean terrain elevation */
   double getBandTerrainHeight(const XMLMetadataSupplier&) const;
 

Modules/Core/Metadata/src/CMakeLists.txt

@@ -77,6 +77,8 @@ set(OTBMetadata_SRC
 
   otbMetadataStorageInterface.cxx
   otbMetadataSupplierInterface.cxx
+
+  otbDateTime.cxx
   )
 
 add_library(OTBMetadata ${OTBMetadata_SRC})
@@ -86,6 +88,7 @@ target_link_libraries(OTBMetadata
   ${OTBOSSIMAdapters_LIBRARIES}
   ${Boost_LIBRARIES}
   ${OTBGdalAdapters_LIBRARIES}
+  ${OTBBoostAdapters_LIBRARIES}
   )
 
 otb_module_target(OTBMetadata)

Modules/Core/Metadata/src/otbCosmoImageMetadataInterface.cxx

@@ -317,6 +317,12 @@ std::vector<Orbit> CosmoImageMetadataInterface::getOrbits(const std::string & re
 
   std::vector<Orbit> orbitVector;
 
+  // Helper function to convert to a two digit string.
+  auto formatTimeInt = [](int i) { return (i < 10 ? "0" + std::to_string(i)
+                                                 : std::to_string(i) );};
+  auto formatTimeDouble = [](double d) { return (d < 10 ? "0" + std::to_string(d)
+                                                       : std::to_string(d) );};
+
   std::ostringstream oss;
   for (int i = 0; i != stateVectorList_size ; ++i)
   {
@@ -328,9 +334,11 @@ std::vector<Orbit> CosmoImageMetadataInterface::getOrbits(const std::string & re
     int hour = static_cast<int> (total_seconds/3600.0);
     int minutes = static_cast<int> ((total_seconds-hour*3600)/60.0);
     double seconds = total_seconds - hour*3600 - minutes*60;
-    std::string timestr = reference_UTC + " " + std::to_string(hour) + ":" + std::to_string(minutes) + ":" + std::to_string(seconds);
+    seconds += 0.5;
+
+    std::string timeStr = reference_UTC + "T" + formatTimeInt(hour) + ":" + formatTimeInt(minutes) + ":" + formatTimeDouble(seconds);
 
-    orbit.time = boost::posix_time::time_from_string(timestr);
+    orbit.time = MetaData::ReadFormattedDate(timeStr);
 
     orbit.position[0] = std::stod(vECEF_sat_pos[i*3 + 0]) ;
     orbit.position[1] = std::stod(vECEF_sat_pos[i*3 + 1]) ;
@@ -356,16 +364,8 @@ std::vector<BurstRecord> CosmoImageMetadataInterface::CreateBurstRecord(const st
   record.endLine = endLine;
   record.endSample = endSample;
 
-
-  std::stringstream ss;
-  auto facet = new boost::posix_time::time_input_facet("%Y-%m-%dT%H:%M:%S%F");
-  ss.imbue(std::locale(std::locale(), facet));
-
-  ss << firstLineTimeStr;
-  ss >> record.azimuthStartTime;
-
-  ss << lastLineTimeStr;
-  ss >> record.azimuthStopTime;
+  record.azimuthStartTime = MetaData::ReadFormattedDate(firstLineTimeStr, "%Y-%m-%dT%H:%M:%S%F");
+  record.azimuthStopTime = MetaData::ReadFormattedDate(lastLineTimeStr, "%Y-%m-%dT%H:%M:%S%F");
 
   record.azimuthAnxTime = 0.;
 
@@ -462,7 +462,7 @@ void CosmoImageMetadataInterface::ParseGdal(ImageMetadata & imd)
   }
   sarParam.rangeSamplingRate = 1 / rangeTimeInterval;
 
-  sarParam.azimuthTimeInterval = boost::posix_time::precise_duration(
+  sarParam.azimuthTimeInterval = MetaData::DurationType(
       std::stod(metadataBands[0]["S01_SBI_Line_Time_Interval"])
       * 1e6);
   sarParam.nearRangeTime = std::stod(metadataBands[0]["S01_SBI_Zero_Doppler_Range_First_Time"]);

Modules/Core/Metadata/src/otbDateTime.cxx

+/*
+ * Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * 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 "otbDateTime.h"
+#include "otbStringUtilities.h"
+
+#if OTB_USE_BOOST_TIME
+namespace boost 
+{ 
+namespace posix_time 
+{
+
+time_duration abs(time_duration d)
+{
+  if(d.is_negative())
+    d = d.invert_sign();
+  return d;
+}
+
+
+}
+}
+
+namespace otb
+{
+namespace MetaData
+{
+
+TimeType ReadFormattedDate(const std::string & dateStr, const std::string & format)
+{
+  MetaData::TimeType outputDate;
+  std::stringstream ss;
+  auto facet = new boost::posix_time::time_input_facet(format);
+  ss.imbue(std::locale(std::locale(), facet));
+  ss << dateStr;
+  ss >> outputDate;
+  return outputDate;
+}
+
+DurationType seconds(double input)
+{
+  return boost::posix_time::precise_duration(input * 1e6);
+}
+
+}
+}
+
+#else //OTB_USE_BOOST_TIME
+
+namespace otb
+{
+namespace MetaData
+{
+
+ModifiedJulianDate ReadFormattedDate(const std::string & dateStr, const std::string & format)
+{
+  otb::string_view fmt  = format;
+  otb::string_view d    = dateStr;
+  int year    = 0;
+  int month   = 0;
+  int day     = 0;
+  int hours   = 0;
+  int minutes = 0;
+  int seconds = 0;
+  double       frac_sec   = 0.0;
+
+  for ( ; !fmt.empty() ; fmt.remove_prefix(1))
+  {
+    switch (fmt.front())
+    {
+      case '%':
+        assert(!fmt.empty()); // invalid format => error
+        fmt.remove_prefix(1);
+        switch (fmt.front()) 
+        {
+          case '%':
+            assert(!d.empty()); // invalid format => error
+            if (fmt.front() != d.front()) 
+            {
+              throw std::runtime_error("Date "+ dateStr +" doesn't match format (" + format + ")");
+            }
+            d.remove_prefix(1);
+            break;
+          case 'Y': year    = otb::details::decode_uint(d); break;
+          case 'm': month   = otb::details::decode_uint(d); break;
+          case 'd': day     = otb::details::decode_uint(d); break;
+          case 'H': hours   = otb::details::decode_uint(d); break;
+          case 'M': minutes = otb::details::decode_uint(d); break;
+          case 'S': seconds = otb::details::decode_uint(d); break;
+          case 'F': // subseconds
+            if (!d.empty() && d.front() == '.')
+            {
+              d.remove_prefix(1);
+              std::size_t nsecs_nb_digits = d.size();
+              frac_sec = otb::details::decode_uint(d);
+              nsecs_nb_digits -= d.size();
+              frac_sec /= std::pow((double)10.0, (int)nsecs_nb_digits);
+            }
+            break;
+          default:
+            throw std::logic_error("Unsupported date format specifier (in "+format+": "+fmt.front()+")");
+         }
+         break;
+      default:
+         if (fmt.front() != d.front()) {
+          throw std::runtime_error("Date "+ dateStr +" doesn't match format (" + format + ")");
+         }
+         d.remove_prefix(1);
+    }
+  }
+
+#if 1
+  // Conversion to julian day
+  // according to http://en.wikipedia.org/wiki/Julian_day
+  // division are integer divisions:
+  int a = (14 - month) / 12;
+  int y = year + 4800 - a;
+  int m = month + 12 * a - 3;
+  int julianDayNumber = day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045;
+
+  // now, the division are NOT integer
+  double julianDay = julianDayNumber +  hours / 24. +  minutes / 1440.+ (seconds + frac_sec) / 86400.;
+
+#else
+  double julianDay = day -32075+1461*(year+4800+(month-14)/12)/4+367*(month-2-(month-14)/12*12)
+           /12-3*((year+4900+(month-14)/12)/100)/4 + (hours/24.0)
+               + (minutes/1440.0)
+               + ((seconds+frac_sec)/86400.0);
+#endif
+
+  // Convert Juian day to modified julian day
+  return ModifiedJulianDate(julianDay - 2400000.5);
+}
+
+}
+}
+
+#endif
\ No newline at end of file

Modules/Core/Metadata/src/otbSarImageMetadataInterface.cxx

@@ -239,10 +239,6 @@ std::vector<Orbit> SarImageMetadataInterface::GetOrbitsGeom() const
   // This streams will hold the iteration number
   std::ostringstream oss;
 
-  std::stringstream ss;
-  auto facet = new boost::posix_time::time_input_facet("%Y-%m-%dT%H:%M:%S%F");
-  ss.imbue(std::locale(std::locale(), facet));
-
   for (int listId = 0 ; listId <= listCount - 1 ; ++listId)
   {
     oss.str("");
@@ -251,8 +247,7 @@ std::vector<Orbit> SarImageMetadataInterface::GetOrbitsGeom() const
     std::string path_root = "orbitList.orbit[" + oss.str() + "]";
     Orbit orbit;
 
-    ss << m_MetadataSupplierInterface->GetAs<std::string>(path_root + ".time");
-    ss >> orbit.time;
+    orbit.time = MetaData::ReadFormattedDate(m_MetadataSupplierInterface->GetAs<std::string>(path_root + ".time"));
 
     orbit.position[0] = m_MetadataSupplierInterface->GetAs<double>(path_root + ".x_pos");
     orbit.position[1] = m_MetadataSupplierInterface->GetAs<double>(path_root + ".y_pos");
@@ -273,10 +268,6 @@ std::vector<BurstRecord> SarImageMetadataInterface::GetBurstRecordsGeom() const
 
   int listCount = m_MetadataSupplierInterface->GetAs<int>(prefix + "geom.bursts.number");
 
-  std::stringstream ss;
-  auto facet = new boost::posix_time::time_input_facet("%Y-%m-%dT%H:%M:%S%F");
-  ss.imbue(std::locale(std::locale(), facet));
-
   const int version = m_MetadataSupplierInterface->GetAs<int>("header.version");