Fixed a bug when the longitude lies in a small sector between the end and start

CMakeLists.txt

@@ -13,6 +13,8 @@ if(${srcdir} STREQUAL ${bindir})
  " Please remove them.")
 endif()
 
+set(STATIC_ANALYZER OFF CACHE BOOL "Using GCC static analyzer if available")
+
 project(odm C CXX)
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
@@ -28,6 +30,12 @@ if(COMPILER_SUPPORTS_DWARF4 AND (CMAKE_BUILD_TYPE STREQUAL "Debug" OR CMAKE_BUIL
  add_compile_options(-gdwarf-4)
 endif()
 
+# GCC static analyzer support check
+CHECK_CXX_COMPILER_FLAG(-fanalyzer COMPILER_SUPPORTS_ANALYZER)
+if(COMPILER_SUPPORTS_ANALYZER AND STATIC_ANALYZER)
+ add_compile_options(-fanalyzer)
+endif()
+
 add_library(teos STATIC GSW-C/gsw_oceanographic_toolbox.c GSW-C/gsw_saar.c)
 set_target_properties(teos PROPERTIES LINKER_LANGUAGE C)
 

include/AVISOLOCAL.h

@@ -0,0 +1,64 @@
+#pragma once
+#include "mregex.h"
+#include "ncfuncs.h"
+#include "simple2ddata.h"
+#include <dirent.h>
+#include <utility>
+
+using michlib::Ceil;
+using michlib::DetGeoDomain;
+using michlib::Floor;
+using michlib::GPL;
+using michlib::int2;
+using michlib::int4;
+using michlib::RegExp;
+
+class AVISOLOCALData: public NCFuncs
+{
+ std::vector<MDateTime> times;
+ MString                datapath;
+
+ struct Parameters: public BaseParameters
+ {
+  real lonb, latb, lone, late;
+  virtual ~Parameters() override = default;
+ };
+
+ static MString Title() { return "AVISO local mirror"; }
+
+ public:
+ using Data = Simple2DData;
+
+ AVISOLOCALData() = default;
+
+ MString Info() const;
+ // TODO: RetVal
+ MString Open(const CLArgs& args);
+
+ std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const;
+
+ bool isOk() const { return times.size() > 0; }
+
+ size_t NTimes() const { return times.size(); }
+
+ MDateTime Time(size_t i) const
+ {
+  if(!isOk() || i >= times.size()) return MDateTime();
+  return times[i];
+ }
+
+ time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
+
+ explicit operator bool() const { return times.size() > 0; }
+
+ bool CheckVar(const MString& vname) const
+ {
+  NCFileA nc;
+  nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc");
+  return NCFuncs::CheckVar(vname, [&nc = std::as_const(nc)](const MString& vn) { return HaveVar(nc, vn); });
+ }
+
+ Data Read(const MString& vname, const BaseParameters* ip, size_t i) const;
+
+ template<class DataType> Data ReadVarRaw(const NCFileA& nc, const MString& name, const struct Parameters* p) const;
+};

include/BINFILE.h

@@ -0,0 +1,44 @@
+#pragma once
+#include "DataAdapters/datafile.cpp"
+#include "DataAdapters/findfilebyid.cpp"
+#include "ParameterList.h"
+#include "mregex.h"
+#include "simple2ddata.h"
+
+using michlib::GPL;
+using michlib::MDateTime;
+using michlib::RegExp;
+
+class BINFILEData
+{
+ std::vector<MDateTime>            times;
+ MString                           datapath;
+ std::unique_ptr<michlib::IntData> data;
+
+ public:
+ using Data = Simple2DData;
+
+ BINFILEData() = default;
+
+ MString Info() const;
+ // TODO: RetVal
+ MString Open(const CLArgs& args);
+
+ bool isOk() const { return times.size() > 0; }
+
+ size_t NTimes() const { return times.size(); }
+
+ MDateTime Time(size_t i) const
+ {
+  if(!isOk() || i >= times.size()) return MDateTime();
+  return times[i];
+ }
+
+ time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
+
+ explicit operator bool() const { return times.size() > 0; }
+
+ bool CheckVar(const MString& vname) const { return vname == "u" || vname == "v" || vname == "U" || vname == "U2"; }
+
+ Data Read(const MString& vname, size_t i) const;
+};

include/ParseArgs.h

@@ -2,6 +2,7 @@
 #include "GPL.h"
 
 using michlib::MString;
+using michlib::SList;
 using CLArgs = std::map<MString, MString>;
 
 CLArgs ParseArgs(int argc, char** argv);

include/actioninfo.h

@@ -0,0 +1,17 @@
+#pragma once
+#include "actions.h"
+
+using michlib::message;
+
+template<class D> struct DoAction_<D, ActionID::INFO>
+{
+ static MString DoAction(const CLArgs& args, D& data);
+};
+
+template<class D> MString DoAction_<D, ActionID::INFO>::DoAction([[maybe_unused]] const CLArgs& args, D& data)
+{
+ auto info = data.Info();
+ if(!info.Exist()) return "No info";
+ message(info);
+ return "";
+};

include/actions.h

@@ -0,0 +1,246 @@
+#pragma once
+#include "basedata.h"
+#include "mdatetime.h"
+#include "mregex.h"
+#include "uvdata.h"
+#include <variant>
+
+using michlib::BFileW;
+using michlib::MDateTime;
+
+using TIndex = std::vector<size_t>;
+
+enum class ActionID
+{
+ INFO,
+ TSC,
+ UV
+};
+
+template<ActionID id> struct ActionCarrier
+{
+ constexpr static ActionID action = id;
+};
+
+template<class T, ActionID id> struct IsActionSupportedTest
+{
+ static constexpr bool ans = false;
+};
+
+template<class T> struct IsActionSupportedTest<T, ActionID::INFO>
+{
+ static constexpr bool ans = InfoSupported<T>;
+};
+template<class T> struct IsActionSupportedTest<T, ActionID::TSC>
+{
+ static constexpr bool ans = ReadPSupported<T> || ReadSupported<T>;
+};
+template<class T> struct IsActionSupportedTest<T, ActionID::UV>
+{
+ static constexpr bool ans = ReadPSupported<T> || ReadSupported<T>;
+};
+
+template<class T, ActionID id> constexpr bool IsActionSupported = IsActionSupportedTest<T, id>::ans;
+
+using ActionVariants = std::variant<ActionCarrier<ActionID::INFO>, ActionCarrier<ActionID::TSC>, ActionCarrier<ActionID::UV>>;
+
+class Action: public ActionVariants
+{
+ public:
+ Action() = default;
+ Action(ActionVariants&& d): ActionVariants(std::move(d)) {}
+ MString Init(const CLArgs& args)
+ {
+  MString act = args.contains("action") ? args.at("action") : "info";
+  if(act == "info")
+   *this = ActionVariants(ActionCarrier<ActionID::INFO>());
+  else if(act == "tsc")
+   *this = ActionVariants(ActionCarrier<ActionID::TSC>());
+  else if(act == "uv")
+   *this = ActionVariants(ActionCarrier<ActionID::UV>());
+  else
+   return "Unknown action: " + act;
+  return "";
+ }
+};
+
+template<class D, ActionID id> struct DoAction_
+{
+ static MString DoAction(const CLArgs& args, D& data) { return "Internal error"; }
+};
+
+template<ActionID id, class D> MString DoAction(const CLArgs& args, D& data) { return DoAction_<D, id>::DoAction(args, data); }
+
+template<class D> size_t GetTIndex(const D& data, const MDateTime& t)
+{
+ size_t nt = data.NTimes();
+
+ if(t <= data.Time(0)) return 0;
+ if(t >= data.Time(nt - 1)) return nt - 1;
+ for(size_t i = 0; i < nt - 1; i++)
+  if(t >= data.Time(i) && t <= data.Time(i + 1)) return (t - data.Time(i) <= data.Time(i + 1) - t) ? i : (i + 1);
+ return 0;
+}
+
+template<class D> std::pair<TIndex, MString> GetTIndexes(const D& data, const CLArgs& args, michlib_internal::ParameterListEx& pars)
+{
+ TIndex tindexes;
+
+ if(args.contains("time") && (args.contains("timeb") || args.contains("timee")))
+  return {tindexes, "Time must be set via time parameter or timeb and timee parameter but not via both"};
+ if(!(args.contains("time") || (args.contains("timeb") && args.contains("timee")))) return {tindexes, "Time must be set via time parameter or timeb and timee parameter"};
+
+ if(args.contains("time"))
+ {
+  MString   regex = args.at("time");
+  MDateTime time;
+
+  if(time.FromString(regex)) return {TIndex(1, GetTIndex(data, time)), ""};             // Time, not regex
+  if(regex == "BEGIN" || regex == "BEG" || regex == "FIRST") return {TIndex(1, 0), ""}; // First time
+  if(regex == "END" || regex == "LAST") return {TIndex(1, data.NTimes() - 1), ""};      // Last time
+
+  michlib::RegExpSimple reg(regex.Buf());
+  if(reg.Compile() != 0) return {tindexes, "Bad regular expression: " + regex};
+
+  for(size_t i = 0; i < data.NTimes(); i++)
+  {
+   MString date = data.Time(i).ToString();
+   if(reg.Match(date.Buf())) tindexes.push_back(i);
+  }
+  if(tindexes.size() == 0) return {tindexes, "There are no times matching the regular expression: " + regex};
+  if(tindexes.size() == 1)
+   pars.SetParameter("time", data.Time(tindexes[0]).ToString());
+  else
+   pars.SetParameter("timeregex", args.at("time"));
+ }
+ else
+ {
+  MString   tb = args.at("timeb");
+  MString   te = args.at("timee");
+  MDateTime b(tb), e(te);
+
+  auto nt = data.NTimes();
+  if(tb == "BEGIN" || tb == "BEG") b = data.Time(0);
+  if(te == "LAST" || te == "END") e = data.Time(nt - 1);
+
+  const MDateTime& beg = (b < e) ? b : e;
+  const MDateTime& end = (b > e) ? b : e;
+
+  if(beg > data.Time(nt - 1)) return {tindexes, "Begin time " + b.ToTString() + " is greater then end time in the dataset " + data.Time(nt - 1).ToTString()};
+  if(end < data.Time(0)) return {tindexes, "End time " + e.ToTString() + " is lesser then begin time in the dataset " + data.Time(0).ToTString()};
+
+  size_t ib = 0, ie = nt - 1;
+  for(size_t i = 0; i < nt; i++)
+   if(data.Time(i) >= beg)
+   {
+    ib = i;
+    break;
+   }
+
+  for(size_t i = nt; i != 0; i--)
+   if(data.Time(i - 1) <= end)
+   {
+    ie = i - 1;
+    break;
+   }
+
+  tindexes.resize(ie - ib + 1);
+  for(size_t i = 0; i < ie - ib + 1; i++) tindexes[i] = i + ib;
+
+  if(tindexes.size() == 0) return {tindexes, "There are no times between " + b.ToString() + " and " + e.ToString()};
+  pars.SetParameter("timeb", b.ToString());
+  pars.SetParameter("timee", e.ToString());
+ }
+
+ return {tindexes, ""};
+}
+
+template<class D> BaseData Read(const D& data, const MString& vname, const BaseParameters* p, const TIndex& tindex)
+{
+ using RT = ReadType<D>;
+ size_t ind;
+ if(tindex.size() == 1)
+ {
+  ind = tindex[0];
+  michlib::message("Time: " + data.Time(ind).ToTString());
+  if constexpr(ReadPSupported<D>)
+   return data.Read(vname, p, ind);
+  else if constexpr(ReadSupported<D>)
+   return data.Read(vname, ind);
+ }
+ else
+ {
+  Averager<RT> out;
+  bool         ok = true;
+  for(size_t i = 0; i < tindex.size(); i++)
+  {
+   if(!ok) break;
+   ind = tindex[i];
+   michlib::message("Time: " + data.Time(ind).ToTString());
+   RT dat;
+   if constexpr(ReadPSupported<D>)
+    dat = data.Read(vname, p, ind);
+   else if constexpr(ReadSupported<D>)
+    dat = data.Read(vname, ind);
+   if(dat)
+    out.Add(dat);
+   else
+    ok = false;
+  }
+  if(ok) return out.Div();
+ }
+ return BaseData();
+}
+
+template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex)
+{
+ using RT = ReadType<D>;
+ using UV = UVData<RT>;
+ size_t ind;
+ if(tindex.size() == 1)
+ {
+  ind = tindex[0];
+  michlib::message("Time: " + data.Time(ind).ToTString());
+  RT u, v;
+  if constexpr(ReadPSupported<D>)
+  {
+   u = data.Read("u", p, ind);
+   v = data.Read("v", p, ind);
+  }
+  else if constexpr(ReadSupported<D>)
+  {
+   u = data.Read("u", ind);
+   v = data.Read("v", ind);
+  }
+  return UV(u, v);
+ }
+ else
+ {
+  Averager<UV> out;
+  bool         ok = true;
+  for(size_t i = 0; i < tindex.size(); i++)
+  {
+   if(!ok) break;
+   ind = tindex[i];
+   michlib::message("Time: " + data.Time(ind).ToTString());
+   RT u, v;
+   if constexpr(ReadPSupported<D>)
+   {
+    u = data.Read("u", p, ind);
+    v = data.Read("v", p, ind);
+   }
+   else if constexpr(ReadSupported<D>)
+   {
+    u = data.Read("u", ind);
+    v = data.Read("v", ind);
+   }
+   UV dat(u, v);
+   if(dat)
+    out.Add(dat);
+   else
+    ok = false;
+  }
+  if(ok) return out.Div();
+ }
+ return UV();
+}

include/actiontsc.h

@@ -0,0 +1,54 @@
+#pragma once
+#include "actions.h"
+
+template<class D> struct DoAction_<D, ActionID::TSC>
+{
+ static MString DoAction(const CLArgs& args, D& data);
+};
+
+template<class D> MString DoAction_<D, ActionID::TSC>::DoAction(const CLArgs& args, D& ds)
+{
+ michlib_internal::ParameterListEx pars;
+ pars.UsePrefix("");
+ pars.SetParameter("source", args.at("source"));
+
+ auto [tindexes, err] = GetTIndexes(ds, args, pars);
+ if(err.Exist()) return err;
+
+ if(!args.contains("var")) return "Variable not specified";
+ MString vname = args.at("var");
+ if(!ds.CheckVar(vname)) return "Variable " + vname + " not exists in this dataset";
+ pars.SetParameter("variable", vname);
+
+ std::unique_ptr<const BaseParameters> sourcepars;
+ if constexpr(ParametersSupported<D>)
+ {
+  auto [p, err] = ds.Parameters(pars, args);
+  if(err.Exist()) return err;
+  sourcepars.reset(p);
+ }
+ auto p = sourcepars.get();
+
+ auto data = Read(ds, vname, p, tindexes);
+ if(!data) return "Can't read data";
+
+ BFileW  fw;
+ MString name = args.contains("out") ? args.at("out") : "";
+ if(!name.Exist()) name = "out.bin";
+
+ fw.Create(name, 3);
+ fw.SetColumnName(1, "Longitude");
+ fw.SetColumnName(2, "Latitude");
+ fw.SetColumnName(3, vname);
+ fw.SetParameters(pars);
+ for(size_t i = 0; i < data.N(); i++)
+ {
+  fw.Write(data.Lon(i));
+  fw.Write(data.Lat(i));
+  fw.Write(data(i) == data.Fillval() ? NAN : data(i));
+ }
+ fw.Finalize();
+ fw.Close();
+
+ return "";
+};

include/actionuv.h

@@ -1,87 +1,30 @@
-#define MICHLIB_NOSOURCE
-#include "GPL.h"
-#include "odm.h"
+#pragma once
+#include "actions.h"
 
-MString Data::ActionUV(const CLArgs& args)
+template<class D> struct DoAction_<D, ActionID::UV>
 {
- if(args.contains("time") && (args.contains("timeb") || args.contains("timee"))) return "Time must be set via time parameter or timeb and timee parameter but not via both";
- if(!(args.contains("time") || (args.contains("timeb") && args.contains("timee")))) return "Time must be set via time parameter or timeb and timee parameter";
+ static MString DoAction(const CLArgs& args, D& data);
+};
 
+template<class D> MString DoAction_<D, ActionID::UV>::DoAction(const CLArgs& args, D& ds)
+{
  michlib_internal::ParameterListEx pars;
  pars.UsePrefix("");
-
  pars.SetParameter("source", args.at("source"));
- TIndex tindexes;
 
- if(args.contains("time"))
- {
-  tindexes = GetTIndexes(args.at("time")); // Regular expression or single date
-  if(tindexes.size() == 0) return "There are no times matching the regular expression";
-  if(tindexes.size() == 1)
-   pars.SetParameter("time", Time(tindexes[0]).ToString());
-  else
-   pars.SetParameter("timeregex", args.at("time"));
- }
- else
- {
-  MDateTime tb(args.at("timeb")), te(args.at("timee"));
-  tindexes = GetTIndexes(tb, te);
-  if(tindexes.size() == 0) return "There are no times between " + tb.ToString() + " and " + te.ToString();
-  pars.SetParameter("timeb", tb.ToString());
-  pars.SetParameter("timee", te.ToString());
- }
+ auto [tindexes, err] = GetTIndexes(ds, args, pars);
+ if(err.Exist()) return err;
 
  std::unique_ptr<const BaseParameters> sourcepars;
+ if constexpr(ParametersSupported<D>)
  {
-  auto ppar = std::visit(
-      [&pars = pars, &args = args](const auto& arg) -> std::pair<const BaseParameters*, MString>
-      {
-       using T = std::decay_t<decltype(arg)>;
-       if constexpr(ParametersSupported<T>)
-        return arg.Parameters(pars, args);
-       else
-        return {nullptr, ""};
-      },
-      *this);
-  if(ppar.second.Exist()) return ppar.second;
-  sourcepars.reset(ppar.first);
+  auto [p, err] = ds.Parameters(pars, args);
+  if(err.Exist()) return err;
+  sourcepars.reset(p);
  }
+ auto p = sourcepars.get();
 
- auto   p = sourcepars.get();
- size_t ind;
- auto   read = [p = p, &ind = ind](const auto& arg) -> UVData
- {
-  using T = std::decay_t<decltype(arg)>;
-  if constexpr(ReadUVPSupported<T>)
-   return arg.ReadUV(p, ind);
-  else
-   return UVData();
- };
-
- UVData data;
- if(tindexes.size() == 1)
- {
-  ind = tindexes[0];
-  michlib::message("Time: " + Time(ind).ToTString());
-  data = std::visit(read, *this);
- }
- else
- {
-  Averager<UVData> out;
-  bool             ok = true;
-  for(size_t i = 0; i < tindexes.size(); i++)
-  {
-   if(!ok) break;
-   ind = tindexes[i];
-   michlib::message("Time: " + Time(ind).ToTString());
-   auto dat = std::visit(read, *this);
-   if(dat)
-    out.Add(dat);
-   else
-    ok = false;
-  }
-  if(ok) data = out.Div();
- }
+ auto data = ReadUV(ds, p, tindexes);
  if(!data) return "Can't read data";
 
  // Main file
@@ -89,7 +32,7 @@ MString Data::ActionUV(const CLArgs& args)
  if(name.Exist())
  {
   BFileW fw;
-  fw.Create(name, 7);
+  fw.Create(name, 9);
   fw.SetColumnName(1, "Longitude");
   fw.SetColumnName(2, "Latitude");
   fw.SetColumnName(3, "u, cm/s");
@@ -97,6 +40,8 @@ MString Data::ActionUV(const CLArgs& args)
   fw.SetColumnName(5, "Div, (cm/s)/km");
   fw.SetColumnName(6, "Rot, (cm/s)/km");
   fw.SetColumnName(7, "Okubo-Weiss parameter, (cm^2/s^2)/km^2");
+  fw.SetColumnName(8, "Kinetic energy, cm^2/s^2");
+  fw.SetColumnName(9, "Eddy kinetic energy, cm^2/s^2");
   fw.SetParameters(pars);
   for(size_t i = 0; i < data.N(); i++)
   {
@@ -107,6 +52,8 @@ MString Data::ActionUV(const CLArgs& args)
    fw.Write(data.Div(i) == data.Fillval() ? NAN : data.Div(i));
    fw.Write(data.Rot(i) == data.Fillval() ? NAN : data.Rot(i));
    fw.Write(data.OW(i) == data.Fillval() ? NAN : data.OW(i));
+   fw.Write(data.U2(i) == data.Fillval() ? NAN : data.U2(i));
+   fw.Write((data.U(i) == data.Fillval() || data.V(i) == data.Fillval()) ? NAN : (data.U2(i) - (data.U(i) * data.U(i) + data.V(i) * data.V(i))));
   }
   fw.Finalize();
   fw.Close();
@@ -166,4 +113,4 @@ MString Data::ActionUV(const CLArgs& args)
  }
 
  return "";
-}
+};

include/basedata.h

@@ -1,9 +1,7 @@
 #pragma once
-#include "comdefs.h"
+#include "traits.h"
 #include <vector>
 
-using michlib::real;
-
 class BaseParameters
 {
  public:
@@ -49,33 +47,17 @@ class BaseData
  explicit operator bool() const { return N() != 0; }
 };
 
-template<class T>
-concept HaveU = requires(T t) {
- {
-  t.template U(0)
- } -> std::convertible_to<real>;
-};
-
-template<class T>
-concept HaveV = requires(T t) {
- {
-  t.template V(0)
- } -> std::convertible_to<real>;
-};
-
-template<class Data> class Averager: public Data
+template<class Data> class DefaultAverager: public Data
 {
+ static constexpr bool isuv = IsUVData<Data>;
  std::vector<size_t>   count;
- static constexpr bool haveu = HaveU<Data>;
- static constexpr bool havev = HaveV<Data>;
 
  static bool DataOk(const Data* d, size_t i)
  {
-  bool uok = true;
-  bool vok = true;
-  if constexpr(haveu) uok = d->U(i) != Data::Fillval();
-  if constexpr(havev) uok = d->V(i) != Data::Fillval();
-  return uok && vok;
+  if constexpr(isuv)
+   return d->U(i) != Data::Fillval() && d->V(i) != Data::Fillval();
+  else
+   return d->V(i) != Data::Fillval();
  }
 
  public:
@@ -96,33 +78,54 @@ template<class Data> class Averager: public Data
    {
     if(DataOk(this, i))
     {
-     if constexpr(haveu) Data::U(i) += d.U(i);
-     if constexpr(havev) Data::V(i) += d.V(i);
+     if constexpr(isuv)
+     {
+      Data::U(i) += d.U(i);
+      Data::U2(i) += d.U2(i);
+     }
+     Data::V(i) += d.V(i);
      count[i]++;
     }
     else
     {
-     if constexpr(haveu) Data::U(i) = d.U(i);
-     if constexpr(havev) Data::V(i) = d.V(i);
-     count[i] = 1;
+     if constexpr(isuv)
+     {
+      Data::U(i)  = d.U(i);
+      Data::U2(i) = d.U2(i);
+     }
+     Data::V(i) = d.V(i);
+     count[i]   = 1;
     }
    }
  }
 
- Averager& Div()
+ DefaultAverager& Div()
  {
   if(!*this) return *this;
   for(size_t i = 0; i < Data::N(); i++)
-  {
-   if constexpr(haveu)
-   {
-    if(count[i] != 0) Data::U(i) /= count[i];
-   }
-   if constexpr(havev)
+   if(count[i] != 0)
    {
-    if(count[i] != 0) Data::V(i) /= count[i];
+    if constexpr(isuv)
+    {
+     Data::U(i) /= count[i];
+     Data::U2(i) /= count[i];
+    }
+    Data::V(i) /= count[i];
    }
-  }
   return *this;
  }
 };
+
+namespace internal
+{
+template<class Data, class Dummy = void> struct GetAverager_
+{
+ using type = DefaultAverager<Data>;
+};
+template<class Data> struct GetAverager_<Data, void_<typename Data::Averager>>
+{
+ using type = typename Data::Averager;
+};
+} // namespace internal
+
+template<class Data> using Averager = internal::GetAverager_<Data>::type;

include/data.h

@@ -0,0 +1,17 @@
+#pragma once
+#include "AVISO.h"
+#include "AVISOLOCAL.h"
+#include "BINFILE.h"
+#include "HYCOM.h"
+#include "NEMO.h"
+#include <variant>
+
+using DataVariants = std::variant<NEMOData, HYCOMData, AVISOData, AVISOLOCALData, BINFILEData>;
+class Data: public DataVariants
+{
+ public:
+ Data() = default;
+ Data(DataVariants&& v): DataVariants(std::move(v)) {}
+
+ MString Init(const CLArgs& args);
+};

include/layereddata.h

@@ -1,25 +1,16 @@
 #pragma once
-#include "DataAdapters/ncfilealt.h"
-#include "ParseArgs.h"
 #include "gsw.h"
-#include "mdatetime.h"
+#include "ncfuncs.h"
 #include "simple2ddata.h"
-#include "uvdata.h"
-#include <algorithm>
 #include <memory>
-#include <set>
 
 using michlib::Ceil;
 using michlib::DetGeoDomain;
 using michlib::Floor;
 using michlib::GPL;
 using michlib::int2;
-using michlib::MDateTime;
-using michlib::MString;
-using michlib::NCFileA;
-using michlib::ToGeoDomain;
 
-class LayeredData
+class LayeredData: public NCFuncs
 {
  public:
  using Data = Simple2DData;
@@ -34,32 +25,22 @@ class LayeredData
   public:
   NC(MString&& newurl): url(std::move(newurl)) { nc.Reset(url + "#cache&noprefetch"); }
 
-  MString ReadTimes()
+  MString ReadTimes(const MString& tname)
   {
    if(!nc) return "File not open";
-   auto time = nc.VR("time");
+   auto time = nc.VR(tname);
    if(!time) return "Can't read times";
 
    MDateTime refdate;
    time_t    step = 0;
    {
-    auto units = nc.Attribute<MString>("time", "units");
+    auto units = nc.Attribute<MString>(tname, "units");
     if(!units) return "Can't read refdate";
-
-    MString rstr;
-    auto    words = michlib::Split_on_words(units);
-    auto    ci    = words.begin();
-    if(ci != words.end())
-    {
-     if(*ci == "hours") step = 3600;
-     if(*ci == "days") step = 3600 * 24;
-     ci++;
-    }
-    if(ci != words.end()) ci++;              // skip "since"
-    if(ci != words.end()) rstr = *ci;        // Day
-    if(ci != words.end()) ci++;
-    if(ci != words.end()) rstr += " " + *ci; // Hours
-    if(!refdate.FromString(rstr)) return "Can't parse " + rstr + " to refdate";
+    auto [rd, st, suc] = Refdate(units);
+    if(!suc) return "Can't parse " + units.Get() + " to refdate";
+    if(st == 0) return "Can't get timestep from string " + units.Get();
+    refdate = rd;
+    step    = st;
    }
 
    times.resize(time.DimLen(0));
@@ -72,6 +53,8 @@ class LayeredData
   MDateTime Begin() const { return times.front(); }
   MDateTime End() const { return times.back(); }
 
+  const NCFileA& Get() const { return nc; }
+
   const std::vector<MDateTime>& Times() const { return times; }
   size_t                        Index(MDateTime tm) const
   {
@@ -94,8 +77,7 @@ class LayeredData
  std::vector<NC>        nc;
  std::vector<real>      depths;
  std::vector<MDateTime> times;
- MString                lonname, latname;
- size_t                 nx, ny;
+ struct CoordNames      dname;
  real                   lonb, latb, lone, late;
  real                   lonstep, latstep;
  MString                title;
@@ -155,8 +137,6 @@ class LayeredData
 
  Data Read(const MString& vname, const BaseParameters* ip, size_t i) const;
 
- UVData ReadUV(const BaseParameters* ip, size_t i) const;
-
  bool isOk() const { return nc.size() > 0; }
 
  explicit operator bool() const { return nc.size() > 0; }
@@ -194,16 +174,7 @@ class LayeredData
 
  bool CheckVar(const MString& vname) const
  {
-  if(!HaveVar(vname))
-  {
-   bool varexist = false;
-   if(vname == "temp" && HaveVar("ptemp") && HaveVar("sal")) varexist = true;
-   if(vname == "ptemp" && HaveVar("temp") && HaveVar("sal")) varexist = true;
-   if(vname == "pdens" && (HaveVar("ptemp") || HaveVar("temp")) && HaveVar("sal")) varexist = true;
-   if((vname == "U" || vname == "U2") && HaveVar("u") && HaveVar("v")) varexist = true;
-   if(!varexist) return false;
-  }
-  return true;
+  return NCFuncs::CheckVar(vname, [this](const MString& vn) { return HaveVar(vn); });
  }
 
  private:
@@ -212,15 +183,7 @@ class LayeredData
  bool HaveVar(const MString& vname) const
  {
   for(size_t i = 0; i < nc.size(); i++)
-  {
-   auto head = nc[i]->Header();
-   for(const auto& v: head.Variables())
-   {
-    auto stname = nc[i]->A<MString>(v.Name(), "standard_name");
-    if(!stname) continue;
-    if(StName2Name(stname) == vname) return true;
-   }
-  }
+   if(NCFuncs::HaveVar(nc[i].Get(), vname)) return true;
   return false;
  }
 
@@ -240,20 +203,4 @@ class LayeredData
   }
   return {"", 0, 0};
  }
-
- static MString StName2Name(const MString& stname)
- {
-  if(stname == "sea_water_potential_temperature") return "ptemp";
-  if(stname == "sea_water_temperature") return "temp";
-  if(stname == "sea_water_salinity") return "sal";
-  if(stname == "ocean_mixed_layer_thickness_defined_by_sigma_theta") return "mld";
-  if(stname == "sea_surface_height_above_geoid") return "ssh";
-  if(stname == "sea_surface_elevation") return "ssh";
-  if(stname == "eastward_sea_water_velocity") return "u";
-  if(stname == "northward_sea_water_velocity") return "v";
-  if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "u";
-  if(stname == "surface_geostrophic_northward_sea_water_velocity") return "v";
-  if(stname == "upward_sea_water_velocity") return "w";
-  return "";
- }
 };

include/ncfuncs.h

@@ -0,0 +1,39 @@
+#pragma once
+#include "DataAdapters/ncfilealt.h"
+#include "mdatetime.h"
+#include <set>
+#include <tuple>
+
+using michlib::MDateTime;
+using michlib::MString;
+using michlib::NCFileA;
+
+class NCFuncs
+{
+ public:
+ struct CoordNames
+ {
+  MString lonname, latname, depthname, timename;
+  size_t  nx, ny, nz, nt;
+ };
+
+ static MString                             StName2Name(const MString& stname);
+ static std::tuple<MDateTime, time_t, bool> Refdate(const MString& refdate);
+ static void                                GetVars(const NCFileA& nc, std::set<MString>& vars);
+ static CoordNames                          GetCNames(const NCFileA& nc);
+ static CoordNames                          GetDNames(const NCFileA& nc);
+ static bool                                HaveVar(const NCFileA& nc, const MString& vname);
+ template<class HV> static bool             CheckVar(const MString& vname, HV hv)
+ {
+  if(!hv(vname))
+  {
+   bool varexist = false;
+   if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = true;
+   if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = true;
+   if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = true;
+   if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = true;
+   if(!varexist) return false;
+  }
+  return true;
+ }
+};

include/odm.h

@@ -1,172 +1,19 @@
 #pragma once
-#include "AVISO.h"
-#include "BFileW.h"
-#include "HYCOM.h"
-#include "NEMO.h"
-#include <utility>
-#include <variant>
+#include "actioninfo.h"
+#include "actiontsc.h"
+#include "actionuv.h"
+#include "data.h"
 
-using michlib::BFileW;
 using michlib::errmessage;
-using michlib::GPL;
-using michlib::message;
 
-using DataVariants = std::variant<NEMOData, HYCOMData, AVISOData>;
-
-template<class T>
-concept InfoSupported = requires(T t) {
- {
-  t.template Info()
- } -> std::convertible_to<MString>;
-};
-
-template<class T>
-concept ParametersSupported = requires(T t, michlib_internal::ParameterListEx& pars, const CLArgs& args) {
- {
-  t.template Parameters(pars, args).first
- } -> std::convertible_to<const BaseParameters*>;
-};
-
-template<class T>
-concept ReadPSupported = requires(T t, const MString& vname, const BaseParameters* ip, size_t i) {
- {
-  t.template Read(vname, ip, i)(0)
- } -> std::convertible_to<real>;
-};
-template<class T>
-concept ReadSupported = requires(T t, const MString& vname, size_t i) {
- {
-  t.template Read(vname, i)(0)
- } -> std::convertible_to<real>;
-};
-
-template<class T>
-concept ReadUVPSupported = requires(T t, const BaseParameters* ip, size_t i) {
- {
-  t.template ReadUV(ip, i).U(0, 0) + t.template ReadUV(ip, i).V(0, 0)
- } -> std::convertible_to<real>;
-};
-
-class Data: public DataVariants
+template<class T, ActionID id> consteval bool MustSupported()
 {
- public:
- using TIndex = std::vector<size_t>;
-
- private:
- TIndex GetTIndexes(const MDateTime& b, const MDateTime& e) const
- {
-  TIndex           out;
-  auto             nt  = NTimes();
-  const MDateTime& beg = (b < e) ? b : e;
-  const MDateTime& end = (b > e) ? b : e;
-
-  if(beg > Time(nt - 1) || end < Time(0)) return out;
-
-  size_t ib = 0, ie = nt - 1;
-  for(size_t i = 0; i < nt; i++)
-   if(Time(i) >= beg)
-   {
-    ib = i;
-    break;
-   }
-
-  for(size_t i = nt; i != 0; i--)
-   if(Time(i - 1) <= end)
-   {
-    ie = i - 1;
-    break;
-   }
-
-  out.resize(ie - ib + 1);
-  for(size_t i = 0; i < ie - ib + 1; i++) out[i] = i + ib;
-  return out;
- }
-
- TIndex GetTIndexes(const MString& regex)
- {
-  {
-   MDateTime time;
-   if(time.FromString(regex)) return TIndex(1, GetTIndex(time));                   // Time, not regex
-   if(regex == "BEGIN" || regex == "BEG" || regex == "FIRST") return TIndex(1, 0); // First time
-   if(regex == "END" || regex == "LAST") return TIndex(1, NTimes() - 1);           // Last time
-  }
-
-  TIndex out;
-
-  std::unique_ptr<regex_t> regbuf(new regex_t);
-  if(0 != regcomp(regbuf.get(), regex.Buf(), REG_EXTENDED | REG_NOSUB)) return out;
-
-  for(size_t i = 0; i < NTimes(); i++)
-  {
-   MString date = Time(i).ToString();
-   if(0 != regexec(regbuf.get(), date.Buf(), 0, 0, 0)) continue;
-   out.push_back(i);
-  }
-
-  return out;
- }
-
- size_t GetTIndex(const MDateTime& t)
- {
-  size_t nt = NTimes();
-
-  if(t <= Time(0)) return 0;
-  if(t >= Time(nt - 1)) return nt - 1;
-  for(size_t i = 0; i < nt - 1; i++)
-   if(t >= Time(i) && t <= Time(i + 1)) return (t - Time(i) <= Time(i + 1) - t) ? i : (i + 1);
-  return 0;
- }
-
- public:
- Data() = default;
-
- Data(DataVariants&& d): DataVariants(std::move(d)) {}
-
- MDateTime Time(size_t it) const
- {
-  return std::visit([it = it](const auto& arg) -> auto { return arg.Time(it); }, *this);
- }
- size_t NTimes() const
- {
-  return std::visit([](const auto& arg) -> auto { return arg.NTimes(); }, *this);
- }
-
- auto CheckVar(const MString& vname) const
- {
-  return std::visit([&vname = std::as_const(vname)](const auto& arg) -> auto { return arg.CheckVar(vname); }, *this);
- }
-
- MString Init(const CLArgs& args)
- {
-  MString src = args.contains("source") ? args.at("source") : "";
-  if(!src.Exist()) return "No source specified";
-  if(src == "NEMO")
-  {
-   NEMOData data;
-   auto     res = data.Open(args);
-   if(res.Exist()) return "Can't open source " + src + ":\n" + res;
-   *this = Data(std::move(data));
-  }
-  else if(src == "AVISO")
-  {
-   AVISOData data;
-   auto      res = data.Open(args);
-   if(res.Exist()) return "Can't open source " + src + ":\n" + res;
-   *this = Data(std::move(data));
-  }
-  else if(src == "HYCOM")
-  {
-   HYCOMData data;
-   auto      res = data.Open(args);
-   if(res.Exist()) return "Can't open source " + src + ":\n" + res;
-   *this = Data(std::move(data));
-  }
-  else
-   return "Unknown source: " + src;
-  return "";
- }
-
- MString ActionInfo(const CLArgs& args);
- MString ActionTsc(const CLArgs& args);
- MString ActionUV(const CLArgs& args);
+ static_assert(IsActionSupported<T, id>);
+ return true;
 };
+
+constexpr bool supported = MustSupported<NEMOData, ActionID::INFO>() && MustSupported<NEMOData, ActionID::TSC>() && MustSupported<NEMOData, ActionID::UV>() &&
+                           MustSupported<HYCOMData, ActionID::INFO>() && MustSupported<HYCOMData, ActionID::TSC>() && MustSupported<HYCOMData, ActionID::UV>() &&
+                           MustSupported<AVISOData, ActionID::INFO>() && MustSupported<AVISOData, ActionID::TSC>() && MustSupported<AVISOData, ActionID::UV>() &&
+                           MustSupported<AVISOLOCALData, ActionID::INFO>() && MustSupported<AVISOLOCALData, ActionID::TSC>() && MustSupported<AVISOLOCALData, ActionID::UV>() &&
+                           MustSupported<BINFILEData, ActionID::INFO>() && MustSupported<BINFILEData, ActionID::TSC>() && MustSupported<BINFILEData, ActionID::UV>();

include/simple2ddata.h

@@ -43,4 +43,7 @@ class Simple2DData: public BaseData
 
  real Lon(size_t i) const { return Lon(i % nx, i / nx); }
  real Lat(size_t i) const { return Lat(i % nx, i / nx); }
+
+ real XStep() const { return xstep; }
+ real YStep() const { return ystep; }
 };

include/traits.h

@@ -0,0 +1,74 @@
+#pragma once
+#include "ParseArgs.h"
+#include <concepts>
+
+class BaseParameters;
+using michlib::real;
+
+template<class T>
+concept InfoSupported = requires(T t) {
+ {
+  t.Info()
+ } -> std::convertible_to<MString>;
+};
+
+template<class T>
+concept ParametersSupported = requires(T t, michlib_internal::ParameterListEx& pars, const CLArgs& args) {
+ {
+  t.Parameters(pars, args).first
+ } -> std::convertible_to<const BaseParameters*>;
+};
+
+template<class T>
+concept ReadPSupported = requires(T t, const MString& vname, const BaseParameters* ip, size_t i) {
+ {
+  t.Read(vname, ip, i)(0)
+ } -> std::convertible_to<real>;
+};
+template<class T>
+concept ReadSupported = requires(T t, const MString& vname, size_t i) {
+ {
+  t.Read(vname, i)(0)
+ } -> std::convertible_to<real>;
+};
+
+template<class T>
+concept ReadUVPSupported = requires(T t, const BaseParameters* ip, size_t i) {
+ {
+  t.ReadUV(ip, i).U(0, 0) + t.template ReadUV(ip, i).V(0, 0)
+ } -> std::convertible_to<real>;
+};
+
+template<class T>
+concept IsUVData = requires(T t) {
+ {
+  t.U(0)
+ } -> std::convertible_to<real>;
+ {
+  t.V(0)
+ } -> std::convertible_to<real>;
+ {
+  t.U2(0)
+ } -> std::convertible_to<real>;
+};
+
+namespace internal
+{
+template<typename...> using void_ = void;
+
+template<class D, bool RP, bool R> struct GetReadType_;
+
+template<class D, bool R> struct GetReadType_<D, true, R>
+{
+ using p    = BaseParameters*;
+ using type = decltype(D().Read(MString(), p(), 0));
+};
+
+template<class D> struct GetReadType_<D, false, true>
+{
+ using type = decltype(D().Read(MString(), 0));
+};
+
+} // namespace internal
+
+template<class D> using ReadType = internal::GetReadType_<D, ReadPSupported<D>, ReadSupported<D>>::type;

include/uvdata.h

@@ -6,22 +6,35 @@ using michlib::M_PI;
 using michlib::real;
 using michlib::Sqrt;
 
+namespace internal
+{
+template<class T>
+concept HaveXYStep = requires(T t) {
+ {
+  t.XStep()
+ } -> std::convertible_to<real>;
+ {
+  t.YStep()
+ } -> std::convertible_to<real>;
+};
+}
+
 enum class Metric
 {
  EUCLID,
  SPHERIC
 };
 
-class UVData
+class UVDataStorage
 {
+ protected:
  static constexpr real fillval = 1.0e10;
- real                  x0 = 0.0, y0 = 0.0;
- size_t                nx = 0, ny = 0;
- real                  xstep = 0.0, ystep = 0.0;
 
- std::vector<real> u, v;
+ std::vector<real> u, v, u2;
+ size_t            nx = 0, ny = 0;
  Metric            metric;
 
+ UVDataStorage() = default;
  real D(real lon1, real lat1, real lon2, real lat2) const
  {
   switch(metric)
@@ -37,43 +50,6 @@ class UVData
  real Ud(size_t ix, size_t iy) const { return IsCoast(ix, iy) ? 0.0 : U(ix, iy); }
  real Vd(size_t ix, size_t iy) const { return IsCoast(ix, iy) ? 0.0 : V(ix, iy); }
 
- real dUdX(size_t ix, size_t iy) const
- {
-  if(IsCoast(ix, iy)) return fillval;
-  if(ix == 0) return (Ud(1, iy) - Ud(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy));
-  if(ix == nx - 1) return (Ud(nx - 1, iy) - Ud(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy));
-  return 0.5 * ((Ud(ix + 1, iy) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
-                (Ud(ix, iy) - Ud(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
- }
- real dVdX(size_t ix, size_t iy) const
- {
-  if(IsCoast(ix, iy)) return fillval;
-  if(ix == 0) return (Vd(1, iy) - Vd(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy));
-  if(ix == nx - 1) return (Vd(nx - 1, iy) - Vd(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy));
-  return 0.5 * ((Vd(ix + 1, iy) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
-                (Vd(ix, iy) - Vd(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
- }
- real dUdY(size_t ix, size_t iy) const
- {
-  if(IsCoast(ix, iy)) return fillval;
-  if(iy == 0) return (Ud(ix, 1) - Ud(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1));
-  if(iy == ny - 1) return (Ud(ix, ny - 1) - Ud(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1));
-  return 0.5 * ((Ud(ix, iy + 1) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
-                (Ud(ix, iy) - Ud(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
- }
- real dVdY(size_t ix, size_t iy) const
- {
-  if(IsCoast(ix, iy)) return fillval;
-  if(iy == 0) return (Vd(ix, 1) - Vd(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1));
-  if(iy == ny - 1) return (Vd(ix, ny - 1) - Vd(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1));
-  return 0.5 * ((Vd(ix, iy + 1) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
-                (Vd(ix, iy) - Vd(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
- }
-
- // For stationary points
- real LonR(real x, [[maybe_unused]] real y) const { return x0 + x * xstep; }
- real LatR([[maybe_unused]] real x, real y) const { return y0 + y * ystep; }
-
  public:
  enum STPOINT
  {
@@ -93,12 +69,6 @@ class UVData
   STPOINT type = NOPOINT;
  };
 
- UVData() = default;
- UVData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_, Metric m_ = Metric::SPHERIC):
-     x0(x0_), y0(y0_), nx(nx_), ny(ny_), xstep(xs_), ystep(ys_), u(nx_ * ny_), v(nx_ * ny_), metric(m_)
- {
- }
-
  const real& U(size_t i) const { return u[i]; }
  const real& V(size_t i) const { return v[i]; }
  const real& U(size_t ix, size_t iy) const { return U(iy * nx + ix); }
@@ -109,38 +79,188 @@ class UVData
  real& U(size_t ix, size_t iy) { return U(iy * nx + ix); }
  real& V(size_t ix, size_t iy) { return V(iy * nx + ix); }
 
+ const real& U2(size_t i) const { return u2[i]; }
+ const real& U2(size_t ix, size_t iy) const { return U2(iy * nx + ix); }
+
+ real& U2(size_t i) { return u2[i]; }
+ real& U2(size_t ix, size_t iy) { return U2(iy * nx + ix); }
+
  size_t N() const { return u.size(); }
  size_t Nx() const { return nx; }
  size_t Ny() const { return ny; }
 
+ explicit operator bool() const { return N() != 0; }
+
+ bool IsCoast(size_t i) const { return U(i) == fillval || V(i) == fillval; }
+ bool IsCoast(size_t ix, size_t iy) const { return U(ix, iy) == fillval || V(ix, iy) == fillval; }
+
+ static real Fillval() { return fillval; }
+};
+
+template<class OneVarData, bool isgrid> class UVDataDims;
+
+template<class OneVarData> class UVDataDims<OneVarData, true>: public UVDataStorage
+{
+ real x0 = 0.0, y0 = 0.0;
+ real xstep = 0.0, ystep = 0.0;
+
+ void Resize(size_t sz)
+ {
+  UVDataStorage::u.resize(sz);
+  UVDataStorage::v.resize(sz);
+  UVDataStorage::u2.resize(sz);
+ }
+
+ protected:
+ UVDataDims() = default;
+ UVDataDims(const OneVarData& u, const OneVarData& v, Metric m = Metric::SPHERIC)
+ {
+  if(!(u && v))
+  {
+   *this = UVDataDims();
+   return;
+  }
+  x0     = u.Lon(0);
+  y0     = u.Lat(0);
+  xstep  = u.XStep();
+  ystep  = u.YStep();
+  metric = m;
+  nx     = u.Nx();
+  ny     = u.Ny();
+
+  Resize(u.N());
+  for(size_t i = 0; i < u.N(); i++)
+  {
+   if(u(i) == u.Fillval() || v(i) == v.Fillval())
+    U(i) = V(i) = U2(i) = Fillval();
+   else
+   {
+    U(i)  = u(i);
+    V(i)  = v(i);
+    U2(i) = u(i) * u(i) + v(i) * v(i);
+   }
+  }
+ }
+
+ public:
  real Lon(size_t ix, [[maybe_unused]] size_t iy) const { return x0 + ix * xstep; }
  real Lat([[maybe_unused]] size_t ix, size_t iy) const { return y0 + iy * ystep; }
-
  real Lon(size_t i) const { return Lon(i % nx, i / nx); }
  real Lat(size_t i) const { return Lat(i % nx, i / nx); }
+};
 
- explicit operator bool() const { return N() != 0; }
+template<class OneVarData> class UVDataDims<OneVarData, false>: public UVDataStorage
+{
+ std::vector<real> lon, lat;
 
- bool IsCoast(size_t i) const { return U(i) == fillval || V(i) == fillval; }
- bool IsCoast(size_t ix, size_t iy) const { return U(ix, iy) == fillval || V(ix, iy) == fillval; }
+ void Resize(size_t sz)
+ {
+  UVDataStorage::u.resize(sz);
+  UVDataStorage::v.resize(sz);
+  UVDataStorage::u2.resize(sz);
+  lon.resize(sz);
+  lat.resize(sz);
+ }
 
- static real Fillval() { return fillval; }
+ protected:
+ UVDataDims() = default;
+ UVDataDims(const OneVarData& u, const OneVarData& v, Metric m = Metric::SPHERIC)
+ {
+  if(!(u && v))
+  {
+   *this = UVDataDims();
+   return;
+  }
+  metric = m;
+  nx     = u.Nx();
+  ny     = u.Ny();
+
+  Resize(u.N());
+  for(size_t i = 0; i < u.N(); i++)
+  {
+   lon[i] = u.Lon(i);
+   lat[i] = u.Lat(i);
+   if(u(i) == u.Fillval() || v(i) == v.Fillval())
+    U(i) = V(i) = U2(i) = Fillval();
+   else
+   {
+    U(i)  = u(i);
+    V(i)  = v(i);
+    U2(i) = u(i) * u(i) + v(i) * v(i);
+   }
+  }
+ }
+
+ public:
+ real Lon(size_t ix, [[maybe_unused]] size_t iy) const { return Lon(nx * iy + ix); }
+ real Lat([[maybe_unused]] size_t ix, size_t iy) const { return Lat(nx * iy + ix); }
+ real Lon(size_t i) const { return lon[i]; }
+ real Lat(size_t i) const { return lat[i]; }
+};
+
+template<class OneVarData> class UVData: public UVDataDims<OneVarData, internal::HaveXYStep<OneVarData>>
+{
+ using B = UVDataDims<OneVarData, internal::HaveXYStep<OneVarData>>;
+ using B::Ud, B::Vd, B::nx, B::ny, B::D;
+
+ // For stationary points. Works only for meridian-parallel grids.
+ real LonR(size_t i, size_t j, real x, [[maybe_unused]] real y) const { return B::Lon(i, j) + x * (B::Lon(i + 1, j) - B::Lon(i, j)); }
+ real LatR(size_t i, size_t j, [[maybe_unused]] real x, real y) const { return B::Lat(i, j) + y * (B::Lat(i, j + 1) - B::Lat(i, j)); }
+
+ real dUdX(size_t ix, size_t iy) const
+ {
+  if(IsCoast(ix, iy)) return Fillval();
+  if(ix == 0) return (Ud(1, iy) - Ud(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy));
+  if(ix == nx - 1) return (Ud(nx - 1, iy) - Ud(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy));
+  return 0.5 * ((Ud(ix + 1, iy) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
+                (Ud(ix, iy) - Ud(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
+ }
+ real dVdX(size_t ix, size_t iy) const
+ {
+  if(IsCoast(ix, iy)) return Fillval();
+  if(ix == 0) return (Vd(1, iy) - Vd(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy));
+  if(ix == nx - 1) return (Vd(nx - 1, iy) - Vd(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy));
+  return 0.5 * ((Vd(ix + 1, iy) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
+                (Vd(ix, iy) - Vd(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
+ }
+ real dUdY(size_t ix, size_t iy) const
+ {
+  if(IsCoast(ix, iy)) return Fillval();
+  if(iy == 0) return (Ud(ix, 1) - Ud(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1));
+  if(iy == ny - 1) return (Ud(ix, ny - 1) - Ud(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1));
+  return 0.5 * ((Ud(ix, iy + 1) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
+                (Ud(ix, iy) - Ud(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
+ }
+ real dVdY(size_t ix, size_t iy) const
+ {
+  if(IsCoast(ix, iy)) return Fillval();
+  if(iy == 0) return (Vd(ix, 1) - Vd(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1));
+  if(iy == ny - 1) return (Vd(ix, ny - 1) - Vd(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1));
+  return 0.5 * ((Vd(ix, iy + 1) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
+                (Vd(ix, iy) - Vd(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
+ }
+
+ public:
+ using B::Fillval, B::IsCoast, B::Lon, B::Lat, B::U, B::V, B::Nx, B::Ny;
+
+ UVData() = default;
+ UVData(const OneVarData& u, const OneVarData& v, Metric m = Metric::SPHERIC): B(u, v, m) {}
 
- real Div(size_t i) const { return Div(i % nx, i / nx); }
- real Rot(size_t i) const { return Rot(i % nx, i / nx); }
+ real Div(size_t i) const { return Div(i % nx, i / B::nx); }
+ real Rot(size_t i) const { return Rot(i % nx, i / B::nx); }
  real Div(size_t ix, size_t iy) const
  {
   if(!*this) return 0.0;
   real ux = dUdX(ix, iy);
   real vy = dVdY(ix, iy);
-  return (ux == fillval || vy == fillval) ? fillval : (ux + vy);
+  return (ux == Fillval() || vy == Fillval()) ? Fillval() : (ux + vy);
  }
  real Rot(size_t ix, size_t iy) const
  {
   if(!*this) return 0.0;
   real vx = dVdX(ix, iy);
   real uy = dUdY(ix, iy);
-  return (vx == fillval || uy == fillval) ? fillval : (vx - uy);
+  return (vx == Fillval() || uy == Fillval()) ? Fillval() : (vx - uy);
  }
  // Okubo-Weiss parameter
  real OW(size_t i) const { return OW(i % nx, i / nx); }
@@ -152,7 +272,7 @@ class UVData
   real vx = dVdX(ix, iy);
   real vy = dVdY(ix, iy);
 
-  if(ux == fillval || uy == fillval || vx == fillval || vy == fillval) return fillval;
+  if(ux == Fillval() || uy == Fillval() || vx == Fillval() || vy == Fillval()) return Fillval();
 
   real sn = ux - vy;
   real ss = vx + uy;
@@ -162,7 +282,7 @@ class UVData
 
  auto StablePoints(size_t ix, size_t iy) const
  {
-  std::vector<struct StPoint> points;
+  std::vector<struct B::StPoint> points;
   if(!*this) return points;
   if(ix >= Nx() - 1 || iy >= Ny() - 1) return points;
   if(IsCoast(ix, iy) || IsCoast(ix + 1, iy) || IsCoast(ix, iy + 1) || IsCoast(ix + 1, iy + 1)) return points;
@@ -202,24 +322,24 @@ class UVData
    if(des > 0.0)
    {
     if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0 && c > 0.0)
-     return UKNOT;
+     return UVDataStorage::UKNOT;
     else if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0)
-     return SADDLE;
+     return UVDataStorage::SADDLE;
     else
-     return SKNOT;
+     return UVDataStorage::SKNOT;
    }
    else
    {
     bool acyclcrit = (av * y + cv > 0.0);
     if(b < 0.0)
-     return acyclcrit ? SACICFOCUS : SCICFOCUS;
+     return acyclcrit ? UVDataStorage::SACICFOCUS : UVDataStorage::SCICFOCUS;
     else
-     return acyclcrit ? UACICFOCUS : UCICFOCUS;
+     return acyclcrit ? UVDataStorage::UACICFOCUS : UVDataStorage::UCICFOCUS;
    }
   };
 
-  if(x1 >= 0.0 && x1 < 1.0 && y1 >= 0.0 && y1 < 1.0) points.emplace_back(LonR(x1 + ix, y1 + iy), LatR(x1 + ix, y1 + iy), PointType(x1, y1));
-  if(x2 >= 0.0 && x2 < 1.0 && y2 >= 0.0 && y2 < 1.0) points.emplace_back(LonR(x2 + ix, y2 + iy), LatR(x2 + ix, y2 + iy), PointType(x2, y2));
+  if(x1 >= 0.0 && x1 < 1.0 && y1 >= 0.0 && y1 < 1.0) points.emplace_back(LonR(ix, iy, x1, y1), LatR(ix, iy, x1, y1), PointType(x1, y1));
+  if(x2 >= 0.0 && x2 < 1.0 && y2 >= 0.0 && y2 < 1.0) points.emplace_back(LonR(ix, iy, x2, y2), LatR(ix, iy, x2, y2), PointType(x2, y2));
   return points;
  }
 };

michlib

@@ -1 +1 @@
-Subproject commit 4294181cfe0dae07f8c46bc3c3c1ea86f377eae8
+Subproject commit a2105413f8a10e8695c85706487d4db5d2f54a61

src/AVISOLOCAL.cpp

@@ -0,0 +1,240 @@
+#define MICHLIB_NOSOURCE
+#include "AVISOLOCAL.h"
+
+MString AVISOLOCALData::Info() const
+{
+ if(!isOk()) return "";
+
+ NCFileA           nc;
+ struct CoordNames cn, dn;
+ std::set<MString> vars;
+
+ nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc");
+ if(!nc) return "Can't open file " + datapath + "/uv-" + times[0].ToString() + ".nc";
+
+ GetVars(nc, vars);
+ dn = GetDNames(nc);
+ cn = GetCNames(nc);
+
+ MString svars;
+ {
+  bool first = true;
+  for(const auto& v: vars)
+  {
+   svars += (first ? "" : ", ") + v;
+   first = false;
+  }
+ }
+
+ auto lons = nc.VR(cn.lonname);
+ auto lats = nc.VR(cn.latname);
+ if(!(lons && lats)) return "Can't get longitudes/latitudes";
+
+ real lonb    = lons(0);
+ real latb    = lats(0);
+ real lone    = lons(dn.nx - 1);
+ real late    = lats(dn.ny - 1);
+ real lonstep = (lone - lonb) / (dn.nx - 1);
+ real latstep = (late - latb) / (dn.ny - 1);
+
+ // clang-format off
+ return
+  "Dataset: " + Title() + "\n" +
+  " Begin date: " + Time(0).ToString() + "\n" +
+  " End date: " + Time(NTimes()-1).ToString() + "\n" +
+  " Time step: " + Timestep() + " seconds\n" +
+  " Time moments: " + NTimes() + "\n" +
+  " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" +
+  " Grid: " + dn.nx + "x" + dn.ny + " (" + lonstep + " x " + latstep + ")\n" +
+  " Supported variables: " + svars;
+ // clang-format on
+}
+
+MString AVISOLOCALData::Open(const CLArgs& args)
+{
+ GPL.UsePrefix("AVISOLOCAL");
+ datapath = GPL.ParameterSValue("Datapath", "");
+ RegExp regex("uv-([0-9]{4}-[0-9]{2}-[0-9]{2}).nc");
+ regex.Compile();
+
+ DIR*           dir = opendir(datapath.Buf());
+ struct dirent* de;
+ if(nullptr == dir) return "Can't open directory " + datapath;
+ while((de = readdir(dir)))
+ {
+  if(!regex.Match(de->d_name)) continue;
+  times.emplace_back(MString(de->d_name + regex.Off(1), regex.Len(1)));
+ }
+ closedir(dir);
+ std::sort(times.begin(), times.end());
+ return "";
+}
+
+std::pair<const BaseParameters*, MString> AVISOLOCALData::Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const
+{
+ std::unique_ptr<struct Parameters> ppar{new struct Parameters};
+ if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"};
+ ppar->lonb = args.at("lonb").ToReal();
+ ppar->lone = args.at("lone").ToReal();
+ ppar->latb = args.at("latb").ToReal();
+ ppar->late = args.at("late").ToReal();
+
+ pars.SetParameter("lonb", ppar->lonb);
+ pars.SetParameter("latb", ppar->latb);
+ pars.SetParameter("lone", ppar->lone);
+ pars.SetParameter("late", ppar->late);
+
+ return {ppar.release(), ""};
+}
+
+AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParameters* ip, size_t i) const
+{
+ if(!isOk()) return Data();
+
+ auto    p = dynamic_cast<const struct Parameters*>(ip);
+ NCFileA nc;
+ MString name    = "";
+ bool    isfloat = false, isint2 = false, isint4 = false;
+
+ nc.Reset(datapath + "/uv-" + times[i].ToString() + ".nc");
+ if(!nc) return Data();
+ {
+  auto head = nc.Header();
+  for(const auto& v: head.Variables())
+  {
+   auto stname = nc.A<MString>(v.Name(), "standard_name");
+   if(!stname) continue;
+   if(StName2Name(stname) == vname)
+   {
+    name    = v.Name();
+    isint2  = v.Type().Id() == NC_SHORT;
+    isint4  = v.Type().Id() == NC_INT;
+    isfloat = v.Type().Id() == NC_FLOAT;
+   }
+  }
+ }
+
+ if(!name.Exist()) // Conversion read
+ {
+  // U and U2 from u and v
+  if(vname == "U" || vname == "U2")
+  {
+   bool square = vname == "U2";
+   auto u      = Read("u", ip, i);
+   auto v      = Read("v", ip, i);
+   if(!(u && v)) return Data();
+   auto out = u;
+   for(size_t ind = 0; ind < out.N(); ind++)
+   {
+    if(u.IsFill(ind) || v.IsFill(ind))
+     out.V(ind) = out.Fillval();
+    else
+     out.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
+   }
+   return out;
+  }
+  return Data();
+ }
+
+ // Direct read
+ if(isint2) return ReadVarRaw<int2>(nc, name, p);
+ if(isint4) return ReadVarRaw<int4>(nc, name, p);
+ if(isfloat) return ReadVarRaw<float>(nc, name, p);
+ return Data();
+}
+
+template<class DataType> AVISOLOCALData::Data AVISOLOCALData::ReadVarRaw(const NCFileA& nc, const MString& name, const struct AVISOLOCALData::Parameters* p) const
+{
+ real     unitmul = 1.0;
+ DataType fill;
+ real     offset = 0.0, scale = 1.0;
+
+ {
+  auto a_fill     = nc.A<DataType>(name, "_FillValue");
+  auto a_offset_d = nc.A<double>(name, "add_offset");
+  auto a_scale_d  = nc.A<double>(name, "scale_factor");
+  auto a_offset_f = nc.A<float>(name, "add_offset");
+  auto a_scale_f  = nc.A<float>(name, "scale_factor");
+  if(!a_fill) return Data();
+  fill = a_fill;
+  if(a_offset_d) offset = a_offset_d;
+  if(a_scale_d) scale = a_scale_d;
+  if(a_offset_f) offset = a_offset_f;
+  if(a_scale_f) scale = a_scale_f;
+ }
+
+ auto cn = GetCNames(nc);
+ auto dn = GetDNames(nc);
+
+ auto lons = nc.VR(cn.lonname);
+ auto lats = nc.VR(cn.latname);
+ if(!(lons && lats)) return Data();
+
+ auto lonb    = lons(0);
+ auto latb    = lats(0);
+ auto lone    = lons(dn.nx - 1);
+ auto late    = lats(dn.ny - 1);
+ auto lonstep = (lone - lonb) / (dn.nx - 1);
+ auto latstep = (late - latb) / (dn.ny - 1);
+
+ auto dom  = DetGeoDomain(lonb, lone);
+ real lon1 = ToGeoDomain(p->lonb, dom);
+ real lon2 = ToGeoDomain(p->lone, dom);
+ real lat1 = p->latb;
+ real lat2 = p->late;
+
+ // Special case when the longitude lies in a small sector between the end and the start
+ if(lon1 < lonb) lon1 = lone;
+ if(lon2 > lone) lon2 = lonb;
+
+ auto yb = static_cast<size_t>(Floor((lat1 - latb) / latstep));
+ auto ye = static_cast<size_t>(Ceil((lat2 - latb) / latstep));
+ if(ye > dn.ny - 1) ye = dn.ny - 1;
+ if(yb >= ye) return Data();
+
+ auto xb = static_cast<size_t>(Floor((lon1 - lonb) / lonstep));
+ auto xe = static_cast<size_t>(Ceil((lon2 - lonb) / lonstep));
+
+ if(xb == xe) return Data();
+
+ auto unit = nc.A<MString>(name, "units");
+ if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0;
+
+ Data data((xb < xe) ? (xe - xb + 1) : (dn.nx + xe - xb + 1), ye - yb + 1, lons(xb), lats(yb), lonstep, latstep);
+
+ if(xb < xe)
+ {
+  auto var = nc.V<DataType>(name, {dn.lonname, xb, xe - xb + 1}, {dn.latname, yb, ye - yb + 1});
+
+  if(!var) return Data();
+  if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data();
+
+  for(size_t ix = 0; ix < var.DimLen(0); ix++)
+   for(size_t iy = 0; iy < var.DimLen(1); iy++)
+   {
+    DataType v   = var(ix, iy);
+    data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+   }
+ }
+ else
+ {
+  auto var1 = nc.V<DataType>(name, {dn.lonname, xb}, {dn.latname, yb, ye - yb + 1});
+  auto var2 = nc.V<DataType>(name, {dn.lonname, 0, xe + 1}, {dn.latname, yb, ye - yb + 1});
+
+  if(!(var1 && var2)) return Data();
+  if((var1.DimLen(0) + var2.DimLen(0)) != data.Nx() || var1.DimLen(1) != data.Ny() || var2.DimLen(1) != data.Ny()) return Data();
+  for(size_t ix = 0; ix < var1.DimLen(0); ix++)
+   for(size_t iy = 0; iy < var1.DimLen(1); iy++)
+   {
+    DataType v   = var1(ix, iy);
+    data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+   }
+  for(size_t ix = 0; ix < var2.DimLen(0); ix++)
+   for(size_t iy = 0; iy < var2.DimLen(1); iy++)
+   {
+    DataType v                    = var2(ix, iy);
+    data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+   }
+ }
+ return data;
+}

src/BINFILE.cpp

@@ -0,0 +1,97 @@
+#define MICHLIB_NOSOURCE
+#include "BINFILE.h"
+
+MString BINFILEData::Info() const
+{
+ if(!isOk()) return "";
+
+ michlib_internal::ParameterListEx pars;
+ pars.AddParameters(data->DataFileParameters());
+ pars.UsePrefix("Datafile_Info");
+
+ MString title  = pars.ParameterSValue("DataSource", "no information");
+ MString lonb   = pars.ParameterSValue("lonb", "");
+ MString lone   = pars.ParameterSValue("lone", "");
+ MString latb   = pars.ParameterSValue("latb", "");
+ MString late   = pars.ParameterSValue("late", "");
+ MString dataid = pars.ParameterSValue("DataID", "no information");
+
+ // clang-format off
+ return
+  "Dataset: " + title + "\n" +
+  " DataId: " + dataid + "\n" +
+  " Begin date: " + Time(0).ToString() + "\n" +
+  " End date: " + Time(NTimes()-1).ToString() + "\n" +
+  " Time step: " + Timestep() + " seconds\n" +
+  " Time moments: " + NTimes() + "\n" +
+  " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" +
+  " Grid: " + data->Nx() + "x" + data->Ny() + "\n" +
+  " Supported variables: u, v, U2";
+ // clang-format on
+}
+
+MString BINFILEData::Open(const CLArgs& args)
+{
+ GPL.UsePrefix("BINFILE");
+ datapath        = GPL.ParameterSValue("Datapath", "");
+ MString dataset = args.contains("dataset") ? args.at("dataset") : "";
+ if(!dataset.Exist()) return "No data file";
+
+ MString file;
+ {
+  RegExp regex("^[0-9a-f]{64}$");
+  regex.Compile();
+  if(regex.Match(dataset.Buf()))
+   file = michlib::FindFileByID(datapath, dataset);
+  else
+   file = (dataset[0] != '/') ? (datapath + "/" + dataset) : dataset;
+ }
+ data.reset(new michlib::IntData);
+ if(!data->Open(file)) return "Can't open file " + file;
+ for(size_t it = 0; it < data->Nt(); it++) times.push_back(data->Date(it));
+
+ return "";
+}
+
+BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const
+{
+ if(!isOk()) return Data();
+
+ // Only rectangular grids are supported
+ real xs = (data->Lon(data->Nx() - 1, data->Ny() - 1) - data->Lon(0, 0)) / (data->Nx() - 1);
+ real ys = (data->Lat(data->Nx() - 1, data->Ny() - 1) - data->Lat(0, 0)) / (data->Ny() - 1);
+ Data out(data->Nx(), data->Ny(), data->Lon(0, 0), data->Lat(0, 0), xs, ys);
+
+ // U and U2 from u and v
+ if(vname == "U" || vname == "U2")
+ {
+  bool square = vname == "U2";
+  auto u      = Read("u", i);
+  auto v      = Read("v", i);
+  if(!(u && v)) return Data();
+  auto out = u;
+  for(size_t ind = 0; ind < out.N(); ind++)
+  {
+   if(u.IsFill(ind) || v.IsFill(ind))
+    out.V(ind) = out.Fillval();
+   else
+    out.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
+  }
+  return out;
+ }
+
+ if(vname == "u" || vname == "v")
+ {
+  bool isu = vname == "u";
+  for(size_t ix = 0; ix < data->Nx(); ix++)
+   for(size_t iy = 0; iy < data->Ny(); iy++)
+   {
+    if(data->IsOk(ix, iy, i))
+     out.V(ix, iy) = isu ? data->Ur(ix, iy, i) : data->Vr(ix, iy, i);
+    else
+     out.V(ix, iy) = out.Fillval();
+   }
+ }
+
+ return out;
+}

src/ParseArgs.cpp

@@ -1,8 +1,6 @@
 #define MICHLIB_NOSOURCE
 #include "ParseArgs.h"
 
-using michlib::SList;
-
 CLArgs ParseArgs(int argc, char** argv)
 {
  CLArgs out;

src/actioninfo.cpp

@@ -1,19 +0,0 @@
-#define MICHLIB_NOSOURCE
-#include "odm.h"
-
-MString Data::ActionInfo(const CLArgs& args)
-{
- MString info = std::visit(
-     [](const auto& arg) -> auto
-     {
-      using T = std::decay_t<decltype(arg)>;
-      if constexpr(InfoSupported<T>)
-       return arg.Info();
-      else
-       return MString();
-     },
-     *this);
- if(!info.Exist()) return "Unsupported combination of action and source";
- message(info);
- return "";
-}

src/actiontsc.cpp

@@ -1,113 +0,0 @@
-#define MICHLIB_NOSOURCE
-#include "GPL.h"
-#include "odm.h"
-
-MString Data::ActionTsc(const CLArgs& args)
-{
- if(args.contains("time") && (args.contains("timeb") || args.contains("timee"))) return "Time must be set via time parameter or timeb and timee parameter but not via both";
- if(!(args.contains("time") || (args.contains("timeb") && args.contains("timee")))) return "Time must be set via time parameter or timeb and timee parameter";
-
- michlib_internal::ParameterListEx pars;
- pars.UsePrefix("");
-
- pars.SetParameter("source", args.at("source"));
- TIndex tindexes;
-
- if(args.contains("time"))
- {
-  tindexes = GetTIndexes(args.at("time")); // Regular expression or single date
-  if(tindexes.size() == 0) return "There are no times matching the regular expression";
-  if(tindexes.size() == 1)
-   pars.SetParameter("time", Time(tindexes[0]).ToString());
-  else
-   pars.SetParameter("timeregex", args.at("time"));
- }
- else
- {
-  MDateTime tb(args.at("timeb")), te(args.at("timee"));
-  tindexes = GetTIndexes(tb, te);
-  if(tindexes.size() == 0) return "There are no times between " + tb.ToString() + " and " + te.ToString();
-  pars.SetParameter("timeb", tb.ToString());
-  pars.SetParameter("timee", te.ToString());
- }
-
- if(!args.contains("var")) return "Variable not specified";
- MString vname = args.at("var");
- if(!CheckVar(vname)) return "Variable " + vname + " not exists in this dataset";
- pars.SetParameter("variable", vname);
-
- std::unique_ptr<const BaseParameters> sourcepars;
- {
-  std::pair<const BaseParameters*, MString> ppar = std::visit(
-      [&pars = pars, &args = args](const auto& arg) -> std::pair<const BaseParameters*, MString>
-      {
-       using T = std::decay_t<decltype(arg)>;
-       if constexpr(ParametersSupported<T>)
-        return arg.Parameters(pars, args);
-       else
-        return {nullptr, ""};
-      },
-      *this);
-  if(ppar.second.Exist()) return ppar.second;
-  sourcepars.reset(ppar.first);
- }
-
- auto   p = sourcepars.get();
- size_t ind;
- auto   read = [&vname = std::as_const(vname), p = p, &ind = ind](const auto& arg) -> BaseData
- {
-  using T = std::decay_t<decltype(arg)>;
-  if constexpr(ParametersSupported<T> && ReadPSupported<T>)
-   return arg.Read(vname, p, ind);
-  else if constexpr(!ParametersSupported<T> && ReadSupported<T>)
-   return arg.Read(vname, ind);
-  else
-   return BaseData();
- };
-
- BaseData data;
- if(tindexes.size() == 1)
- {
-  ind = tindexes[0];
-  michlib::message("Time: " + Time(ind).ToTString());
-  data = std::visit(read, *this);
- }
- else
- {
-  Averager<BaseData> out;
-  bool               ok = true;
-  for(size_t i = 0; i < tindexes.size(); i++)
-  {
-   if(!ok) break;
-   ind = tindexes[i];
-   michlib::message("Time: " + Time(ind).ToTString());
-   auto dat = std::visit(read, *this);
-   if(dat)
-    out.Add(dat);
-   else
-    ok = false;
-  }
-  if(ok) data = out.Div();
- }
- if(!data) return "Can't read data";
-
- BFileW  fw;
- MString name = args.contains("out") ? args.at("out") : "";
- if(!name.Exist()) name = "out.bin";
-
- fw.Create(name, 3);
- fw.SetColumnName(1, "Longitude");
- fw.SetColumnName(2, "Latitude");
- fw.SetColumnName(3, vname);
- fw.SetParameters(pars);
- for(size_t i = 0; i < data.N(); i++)
- {
-  fw.Write(data.Lon(i));
-  fw.Write(data.Lat(i));
-  fw.Write(data(i) == data.Fillval() ? NAN : data(i));
- }
- fw.Finalize();
- fw.Close();
-
- return "";
-}

src/basedata.cpp

@@ -1,3 +1,4 @@
+#define MICHLIB_NOSOURCE
 #include "basedata.h"
 
 BaseParameters::~BaseParameters() {}

src/data.cpp

@@ -0,0 +1,46 @@
+#define MICHLIB_NOSOURCE
+#include "data.h"
+
+MString Data::Init(const CLArgs& args)
+{
+ MString src = args.contains("source") ? args.at("source") : "";
+ if(!src.Exist()) return "No source specified";
+ if(src == "NEMO")
+ {
+  NEMOData data;
+  auto     res = data.Open(args);
+  if(res.Exist()) return "Can't open source " + src + ":\n" + res;
+  *this = Data(std::move(data));
+ }
+ else if(src == "AVISO")
+ {
+  AVISOData data;
+  auto      res = data.Open(args);
+  if(res.Exist()) return "Can't open source " + src + ":\n" + res;
+  *this = Data(std::move(data));
+ }
+ else if(src == "AVISOLOCAL")
+ {
+  AVISOLOCALData data;
+  auto           res = data.Open(args);
+  if(res.Exist()) return "Can't open source " + src + ":\n" + res;
+  *this = Data(std::move(data));
+ }
+ else if(src == "HYCOM")
+ {
+  HYCOMData data;
+  auto      res = data.Open(args);
+  if(res.Exist()) return "Can't open source " + src + ":\n" + res;
+  *this = Data(std::move(data));
+ }
+ else if(src == "BINFILE")
+ {
+  BINFILEData data;
+  auto      res = data.Open(args);
+  if(res.Exist()) return "Can't open source " + src + ":\n" + res;
+  *this = Data(std::move(data));
+ }
+ else
+  return "Unknown source: " + src;
+ return "";
+}

src/layereddata.cpp

@@ -8,22 +8,7 @@ MString LayeredData::Info() const
  for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")";
 
  std::set<MString> vars;
- for(const auto& f: nc)
- {
-  auto head = f->Header();
-  for(const auto& v: head.Variables())
-  {
-   auto ret = f->A<MString>(v.Name(), "standard_name");
-   if(!ret) continue;
-   if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret));
-  }
- }
- if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens");
- if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp");
- if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp");
-
- if(vars.contains("u") && vars.contains("v")) vars.emplace("U");
- if(vars.contains("u") && vars.contains("v")) vars.emplace("U2");
+ for(const auto& f: nc) GetVars(f.Get(), vars);
 
  MString svars;
  {
@@ -43,7 +28,7 @@ MString LayeredData::Info() const
   " Time step: " + Timestep() + " seconds\n" +
   " Time moments: " + NTimes() + "\n" +
   " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" +
-  " Grid: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" +
+  " Grid: " + dname.nx + "x" + dname.ny + " (" + lonstep + " x " + latstep + ")\n" +
   " Depths:" + d + "\n" +
   " Supported variables: " + svars;
  // clang-format on
@@ -69,12 +54,6 @@ MString LayeredData::Open(const MString& dataset)
     nc.clear();
     return "Can't connect to url " + url;
    }
-   MString ret = nc.back().ReadTimes();
-   if(ret.Exist())
-   {
-    nc.clear();
-    return "Can't connect to url " + url + ": " + ret;
-   }
   }
   else
    break;
@@ -82,60 +61,66 @@ MString LayeredData::Open(const MString& dataset)
  }
  if(nc.size() == 0) return "No urls for dataset " + dataset + " specified in config";
 
- auto head = nc[0]->Header();
- lonname = latname = "";
- for(const auto& dim: head.Dimensions())
+ dname = GetDNames(nc[0].Get());
+ if(!(dname.lonname.Exist() && dname.latname.Exist()))
  {
-  if(dim.Name() == "lon" || dim.Name() == "longitude")
-  {
-   lonname = dim.Name();
-   nx      = dim.Len();
-  }
-  if(dim.Name() == "lat" || dim.Name() == "latitude")
-  {
-   latname = dim.Name();
-   ny      = dim.Len();
-  }
+  nc.clear();
+  return "Can't find longitude/latitude";
  }
- if(!(lonname.Exist() && latname.Exist()))
+ if(!dname.timename.Exist())
  {
   nc.clear();
-  return "Can't find longitude/latitude";
+  return "Can't find time";
  }
 
+ auto cn = GetCNames(nc[0].Get());
+
  // Read times
+ for(auto& f: nc)
  {
-  for(const auto& f: nc) times.insert(times.end(), f.Times().begin(), f.Times().end());
-  std::sort(times.begin(), times.end());
-  auto last = std::unique(times.begin(), times.end());
-  times.erase(last, times.end());
+  MString ret = f.ReadTimes(cn.timename);
+  if(ret.Exist())
+  {
+   nc.clear();
+   return ret;
+  }
+  times.insert(times.end(), f.Times().begin(), f.Times().end());
  }
+ std::sort(times.begin(), times.end());
+ auto last = std::unique(times.begin(), times.end());
+ times.erase(last, times.end());
 
- auto rdepths = nc[0]->VR("depth");
- if(rdepths)
+ if(cn.depthname.Exist())
  {
+  auto rdepths = nc[0]->VR(cn.depthname);
+  if(!rdepths)
+  {
+   nc.clear();
+   return "Can't read depths";
+  }
   depths.resize(rdepths.DimLen(0));
   for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i);
  }
  else // Surface only data
  {
   depths.resize(1);
-  depths[0]=0;
+  depths[0] = 0;
  }
 
- auto lons = nc[0]->VR(lonname);
- auto lats = nc[0]->VR(latname);
+ auto lons = nc[0]->VR(cn.lonname);
+ auto lats = nc[0]->VR(cn.latname);
  if(!(lons && lats))
  {
   nc.clear();
   return "Can't get longitudes/latitudes";
  }
+
  lonb    = lons(0);
  latb    = lats(0);
- lone    = lons(nx - 1);
- late    = lats(ny - 1);
- lonstep = (lone - lonb) / (nx - 1);
- latstep = (late - latb) / (ny - 1);
+ lone    = lons(dname.nx - 1);
+ late    = lats(dname.ny - 1);
+ lonstep = (lone - lonb) / (dname.nx - 1);
+ latstep = (late - latb) / (dname.ny - 1);
 
  return "";
 }
@@ -144,7 +129,7 @@ std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_intern
 {
  std::unique_ptr<struct Parameters> ppar{new struct Parameters};
 
- if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger<size_t>();
+ ppar->layer = args.contains("layer") ? args.at("layer").ToInteger<size_t>() : 0;
  if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"};
  real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer);
 
@@ -157,9 +142,13 @@ std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_intern
   real lat1 = args.at("latb").ToReal();
   real lat2 = args.at("late").ToReal();
 
+  // Special case when the longitude lies in a small sector between the end and the start
+  if(lon1 < lonb) lon1 = lone;
+  if(lon2 > lone) lon2 = lonb;
+
   ppar->yb = static_cast<size_t>(Floor((lat1 - latb) / latstep));
   ppar->ye = static_cast<size_t>(Ceil((lat2 - latb) / latstep));
-  if(ppar->ye > ny - 1) ppar->ye = ny - 1;
+  if(ppar->ye > dname.ny - 1) ppar->ye = dname.ny - 1;
   if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"};
 
   ppar->xb = static_cast<size_t>(Floor((lon1 - lonb) / lonstep));
@@ -277,34 +266,14 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters*
   {
    if(v.Dimensions().size() == 3) nodepth = true;
    if(v.Type().Id() == NC_SHORT) return ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
+   if(v.Type().Id() == NC_INT) return ReadVarRaw<int>(nc[id], name, tid, nodepth, p);
    if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(nc[id], name, tid, nodepth, p);
+   if(v.Type().Id() == NC_DOUBLE) return ReadVarRaw<double>(nc[id], name, tid, nodepth, p);
   }
 
  return Data();
 }
 
-UVData LayeredData::ReadUV(const BaseParameters* ip, size_t i) const
-{
- if(!isOk()) return UVData();
-
- auto u = Read("u", ip, i);
- auto v = Read("v", ip, i);
- if(!(u && v)) return UVData();
-
- UVData out{u.Nx(), u.Ny(), u.Lon(0, 0), u.Lat(0, 0), lonstep, latstep};
- for(size_t i = 0; i < out.N(); i++)
- {
-  if(u(i) == Data::Fillval() || v(i) == Data::Fillval())
-   out.U(i) = out.V(i) = UVData::Fillval();
-  else
-  {
-   out.U(i) = u(i);
-   out.V(i) = v(i);
-  }
- }
- return out;
-}
-
 template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct LayeredData::Parameters* p) const
 {
  real     unitmul = 1.0;
@@ -317,8 +286,15 @@ template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f,
   auto a_scale_d  = f->A<double>(name, "scale_factor");
   auto a_offset_f = f->A<float>(name, "add_offset");
   auto a_scale_f  = f->A<float>(name, "scale_factor");
-  if(!a_fill) return Data();
-  fill = a_fill;
+  if(a_fill)
+   fill = a_fill;
+  else
+  {
+   if constexpr(std::is_floating_point_v<DataType>)
+    fill = NAN;
+   else
+    fill = -1;
+  }
   if(a_offset_d) offset = a_offset_d;
   if(a_scale_d) scale = a_scale_d;
   if(a_offset_f) offset = a_offset_f;
@@ -328,12 +304,13 @@ template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f,
  auto unit = f->A<MString>(name, "units");
  if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0;
 
- Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep);
+ Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (dname.nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep);
 
  if(p->xb < p->xe)
  {
-  auto var = nodepth ? f->V<DataType>(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
-                     : f->V<DataType>(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
+  auto var =
+      nodepth ? f->V<DataType>(name, {dname.lonname, p->xb, p->xe - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
+              : f->V<DataType>(name, {dname.lonname, p->xb, p->xe - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, p->layer, 1});
   if(!var) return Data();
   if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data();
 
@@ -341,28 +318,28 @@ template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f,
    for(size_t iy = 0; iy < var.DimLen(1); iy++)
    {
     DataType v   = var(ix, iy);
-    data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+    data(ix, iy) = (v == fill || isnan(v)) ? Data::Fillval() : ((v * scale + offset) * unitmul);
    }
  }
  else
  {
-  auto var1 = nodepth ? f->V<DataType>(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
-                      : f->V<DataType>(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
-  auto var2 = nodepth ? f->V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
-                      : f->V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
+  auto var1 = nodepth ? f->V<DataType>(name, {dname.lonname, p->xb}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
+                      : f->V<DataType>(name, {dname.lonname, p->xb}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, p->layer, 1});
+  auto var2 = nodepth ? f->V<DataType>(name, {dname.lonname, 0, p->xe + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
+                      : f->V<DataType>(name, {dname.lonname, 0, p->xe + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, p->layer, 1});
   if(!(var1 && var2)) return Data();
   if((var1.DimLen(0) + var2.DimLen(0)) != data.Nx() || var1.DimLen(1) != data.Ny() || var2.DimLen(1) != data.Ny()) return Data();
   for(size_t ix = 0; ix < var1.DimLen(0); ix++)
    for(size_t iy = 0; iy < var1.DimLen(1); iy++)
    {
     DataType v   = var1(ix, iy);
-    data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+    data(ix, iy) = (v == fill || isnan(v)) ? Data::Fillval() : ((v * scale + offset) * unitmul);
    }
   for(size_t ix = 0; ix < var2.DimLen(0); ix++)
    for(size_t iy = 0; iy < var2.DimLen(1); iy++)
    {
     DataType v                    = var2(ix, iy);
-    data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+    data(ix + var1.DimLen(0), iy) = (v == fill || isnan(v)) ? Data::Fillval() : ((v * scale + offset) * unitmul);
    }
  }
  return data;

src/ncfuncs.cpp

@@ -0,0 +1,148 @@
+#define MICHLIB_NOSOURCE
+#include "ncfuncs.h"
+
+NCFuncs::CoordNames NCFuncs::GetDNames(const NCFileA& nc)
+{
+ CoordNames out;
+ auto       head = nc.Header();
+ for(const auto& dim: head.Dimensions())
+ {
+  if(dim.Name() == "lon" || dim.Name() == "longitude")
+  {
+   out.lonname = dim.Name();
+   out.nx      = dim.Len();
+  }
+  if(dim.Name() == "lat" || dim.Name() == "latitude")
+  {
+   out.latname = dim.Name();
+   out.ny      = dim.Len();
+  }
+  if(dim.Name() == "depth")
+  {
+   out.depthname = dim.Name();
+   out.nz        = dim.Len();
+  }
+  if(dim.Name() == "time")
+  {
+   out.timename = dim.Name();
+   out.nt       = dim.Len();
+  }
+ }
+ return out;
+}
+
+NCFuncs::CoordNames NCFuncs::GetCNames(const NCFileA& nc)
+{
+ CoordNames out;
+ auto       head = nc.Header();
+ for(const auto& v: head.Variables()) // Try to define coordinates by attribute standard_name or attribute axis
+ {
+  auto stname = nc.A<MString>(v.Name(), "standard_name");
+  auto axis   = nc.A<MString>(v.Name(), "axis");
+  bool islon = false, islat = false, isdepth = false, istime = false;
+  if(!(stname || axis)) continue;
+  if(stname && stname.Get() == "longitude") islon = true;
+  if(stname && stname.Get() == "latitude") islat = true;
+  if(stname && stname.Get() == "depth") isdepth = true;
+  if(stname && stname.Get() == "time") istime = true;
+
+  if(!out.lonname.Exist() && axis && axis.Get() == "X") islon = true;
+  if(!out.latname.Exist() && axis && axis.Get() == "Y") islat = true;
+  if(!out.depthname.Exist() && axis && axis.Get() == "Z") isdepth = true;
+  if(!out.timename.Exist() && axis && axis.Get() == "T") istime = true;
+
+  if(islon) out.lonname = v.Name();
+  if(islat) out.latname = v.Name();
+  if(isdepth) out.depthname = v.Name();
+  if(istime) out.timename = v.Name();
+
+  if(islon) out.nx = v.Dimensions().size();
+  if(islat) out.ny = v.Dimensions().size();
+  if(isdepth) out.nz = v.Dimensions().size();
+  if(istime) out.nt = v.Dimensions().size();
+ }
+ return out;
+}
+
+void NCFuncs::GetVars(const NCFileA& nc, std::set<MString>& vars)
+{
+ auto head = nc.Header();
+ for(const auto& v: head.Variables())
+ {
+  auto ret = nc.A<MString>(v.Name(), "standard_name");
+  if(!ret) continue;
+  if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret));
+ }
+ if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens");
+ if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp");
+ if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp");
+
+ if(vars.contains("u") && vars.contains("v")) vars.emplace("U");
+ if(vars.contains("u") && vars.contains("v")) vars.emplace("U2");
+}
+
+std::tuple<MDateTime, time_t, bool> NCFuncs::Refdate(const MString& refdate)
+{
+ MDateTime out;
+ time_t    step = 0;
+
+ MString rstr;
+ auto    words = michlib::Split_on_words(refdate);
+ auto    ci    = words.begin();
+ if(ci != words.end())
+ {
+  if(*ci == "hours") step = 3600;
+  if(*ci == "days") step = 3600 * 24;
+  ci++;
+ }
+ if(ci != words.end()) ci++;       // skip "since"
+ if(ci != words.end()) rstr = *ci; // Day
+ if(ci != words.end()) ci++;
+ if(ci != words.end()) rstr += " " + *ci; // Hours
+ bool success = out.FromString(rstr);
+ return {out, step, success};
+}
+
+MString NCFuncs::StName2Name(const MString& stname)
+{
+ if(stname == "sea_water_potential_temperature") return "ptemp";
+ if(stname == "sea_water_temperature") return "temp";
+ if(stname == "sea_water_salinity") return "sal";
+ if(stname == "ocean_mixed_layer_thickness_defined_by_sigma_theta") return "mld";
+ if(stname == "sea_surface_height_above_geoid") return "ssh";
+ if(stname == "sea_surface_elevation") return "ssh";
+ if(stname == "eastward_sea_water_velocity") return "u";
+ if(stname == "northward_sea_water_velocity") return "v";
+ if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "u";
+ if(stname == "surface_geostrophic_northward_sea_water_velocity") return "v";
+ if(stname == "upward_sea_water_velocity") return "w";
+ return "";
+}
+
+bool NCFuncs::HaveVar(const NCFileA& nc, const MString& vname)
+{
+ auto head = nc.Header();
+ for(const auto& v: head.Variables())
+ {
+  auto stname = nc.A<MString>(v.Name(), "standard_name");
+  if(!stname) continue;
+  if(StName2Name(stname) == vname) return true;
+ }
+ return false;
+}
+
+/*
+template<class HV> bool NCFuncs::CheckVar(const MString& vname, HV hv)
+{
+ if(!hv(vname))
+ {
+  bool varexist = false;
+  if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = true;
+  if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = true;
+  if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = true;
+  if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = true;
+  if(!varexist) return false;
+ }
+ return true;
+}
+*/

src/odm.cpp

@@ -4,9 +4,9 @@ inline void Usage(const MString& arg0)
 {
  message(arg0 + " (<key>=<value>...)");
  message("Keys are:");
- message(" action. What the program should do. May be: info, tsc, uv. Default: tsc.");
+ message(" action. What the program should do. May be: info, tsc, uv. Default: info.");
  message(" Keys for action=info. Print some information about dataset.");
- message("  source. Required. May be: NEMO, HYCOM, AVISO");
+ message("  source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL, BINFILE");
  message("  Keys for source=NEMO");
  message("   dataset. Can be DT, NRT or NRT6. Default: DT");
  message("  Keys for source=HYCOM");
@@ -14,7 +14,7 @@ inline void Usage(const MString& arg0)
  message("  Keys for source=AVISO");
  message("   dataset. Can be DT, NRT, EckmanDT or EckmanNRT. Default: DT");
  message(" Keys for action=tsc. Get temperature, salinity, chlorofill from dataset.");
- message("  source. Required. May be: NEMO, HYCOM, AVISO");
+ message("  source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL, BINFILE");
  message("  var. Required. May be: U, U2, u, v, temp, ptemp, pdens, sal, chl, mld, ssh or w.");
  message("  time. Time moment or regular expression. If present, timeb and timee must be absent");
  message("  timeb, timee. Time interval. If present, time must be absent");
@@ -29,8 +29,10 @@ inline void Usage(const MString& arg0)
  message("  Keys for source=AVISO");
  message("   dataset. Can be DT, NRT, EckmanDT or EckmanNRT. Default: DT");
  message("   layer and/or depth. Layer or depth of AVISO dataset. If depth is specified, layer is ignored. Both ignored for datasets DT and NRT. Default: layer=0");
+ message("  Keys for source=BINFILE");
+ message("   dataset. Path or DataID of interpolation file");
  message(" Keys for action=uv. Get velocity field and its derivatives.");
- message("  source. Required. May be: NEMO, HYCOM, AVISO");
+ message("  source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL, BINFILE");
  message("  time. Time moment or regular expression. If present, timeb and timee must be absent");
  message("  timeb, timee. Time interval. If present, time must be absent");
  message("  out. Output file for components of velocity field, divergency, rotor and Okubo-Weiss parameter. If absent, this data not calculated.");
@@ -50,7 +52,7 @@ int main(int argc, char** argv)
  }
  auto args = ParseArgs(argc, argv);
 
- MString action = args.contains("action") ? args["action"] : "tsc";
+ MString action = args.contains("action") ? args["action"] : "info";
 
  Data data;
  {
@@ -61,22 +63,32 @@ int main(int argc, char** argv)
    return 1;
   }
  }
-
- MString ret;
- if(action == "info")
-  ret = data.ActionInfo(args);
- else if(action == "tsc")
-  ret = data.ActionTsc(args);
- else if(action == "uv")
- {
-  args["var"] = "U";
-  ret         = data.ActionUV(args);
- }
- else
+ Action act;
  {
-  errmessage("Unknown action " + action);
-  return 1;
+  auto ret = act.Init(args);
+  if(ret.Exist())
+  {
+   errmessage(ret);
+   return 1;
+  }
  }
+
+ auto ret = std::visit(
+     [&act = std::as_const(act), &args = std::as_const(args)](auto& arg)
+     {
+      using DT = std::decay_t<decltype(arg)>;
+      return std::visit(
+          [&data = arg, &args = std::as_const(args)](const auto& arg)
+          {
+           using AT = std::decay_t<decltype(arg)>;
+           if constexpr(IsActionSupported<DT, AT::action>)
+            return DoAction<AT::action>(args, data);
+           else
+            return MString("Unsupported combination of action and source");
+          },
+          act);
+     },
+     data);
  if(ret.Exist())
  {
   errmessage(ret);