Browse Source

Rewrite abstractions

interpolate
Michael Uleysky 2 years ago
parent
commit
f9da82612d
  1. 7
      include/AVISOLOCAL.h
  2. 1
      include/ParseArgs.h
  3. 17
      include/actioninfo.h
  4. 246
      include/actions.h
  5. 54
      include/actiontsc.h
  6. 95
      include/actionuv.h
  7. 81
      include/basedata.h
  8. 16
      include/data.h
  9. 13
      include/layereddata.h
  10. 1
      include/ncfuncs.h
  11. 186
      include/odm.h
  12. 74
      include/traits.h
  13. 252
      include/uvdata.h
  14. 22
      src/AVISOLOCAL.cpp
  15. 2
      src/ParseArgs.cpp
  16. 19
      src/actioninfo.cpp
  17. 113
      src/actiontsc.cpp
  18. 1
      src/basedata.cpp
  19. 39
      src/data.cpp
  20. 22
      src/layereddata.cpp
  21. 4
      src/ncfuncs.cpp
  22. 39
      src/odm.cpp

7
include/AVISOLOCAL.h

@ -1,10 +1,7 @@
#pragma once #pragma once
#include "GPL.h"
#include "ParseArgs.h"
#include "mregex.h" #include "mregex.h"
#include "ncfuncs.h" #include "ncfuncs.h"
#include "simple2ddata.h" #include "simple2ddata.h"
#include "uvdata.h"
#include <dirent.h> #include <dirent.h>
#include <utility> #include <utility>
@ -14,8 +11,6 @@ using michlib::Floor;
using michlib::GPL; using michlib::GPL;
using michlib::int2; using michlib::int2;
using michlib::int4; using michlib::int4;
using michlib::MString;
using michlib::real;
using michlib::RegExp; using michlib::RegExp;
class AVISOLOCALData: public NCFuncs class AVISOLOCALData: public NCFuncs
@ -66,6 +61,4 @@ class AVISOLOCALData: public NCFuncs
Data Read(const MString& vname, const BaseParameters* ip, size_t i) const; 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; template<class DataType> Data ReadVarRaw(const NCFileA& nc, const MString& name, const struct Parameters* p) const;
UVData ReadUV(const BaseParameters* ip, size_t i) const;
}; };

1
include/ParseArgs.h

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

17
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 "";
};

246
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();
}

54
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 "";
};

95
src/actionuv.cpp → include/actionuv.h

@ -1,87 +1,30 @@
#define MICHLIB_NOSOURCE #pragma once
#include "GPL.h" #include "actions.h"
#include "odm.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"; static MString DoAction(const CLArgs& args, D& data);
if(!(args.contains("time") || (args.contains("timeb") && args.contains("timee")))) return "Time must be set via time parameter or timeb and timee parameter"; };
template<class D> MString DoAction_<D, ActionID::UV>::DoAction(const CLArgs& args, D& ds)
{
michlib_internal::ParameterListEx pars; michlib_internal::ParameterListEx pars;
pars.UsePrefix(""); pars.UsePrefix("");
pars.SetParameter("source", args.at("source")); pars.SetParameter("source", args.at("source"));
TIndex tindexes;
if(args.contains("time")) auto [tindexes, err] = GetTIndexes(ds, args, pars);
{ if(err.Exist()) return err;
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());
}
std::unique_ptr<const BaseParameters> sourcepars; std::unique_ptr<const BaseParameters> sourcepars;
if constexpr(ParametersSupported<D>)
{ {
auto ppar = std::visit( auto [p, err] = ds.Parameters(pars, args);
[&pars = pars, &args = args](const auto& arg) -> std::pair<const BaseParameters*, MString> if(err.Exist()) return err;
{ sourcepars.reset(p);
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();
auto p = sourcepars.get(); auto data = ReadUV(ds, p, tindexes);
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();
}
if(!data) return "Can't read data"; if(!data) return "Can't read data";
// Main file // Main file
@ -89,7 +32,7 @@ MString Data::ActionUV(const CLArgs& args)
if(name.Exist()) if(name.Exist())
{ {
BFileW fw; BFileW fw;
fw.Create(name, 7); fw.Create(name, 9);
fw.SetColumnName(1, "Longitude"); fw.SetColumnName(1, "Longitude");
fw.SetColumnName(2, "Latitude"); fw.SetColumnName(2, "Latitude");
fw.SetColumnName(3, "u, cm/s"); 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(5, "Div, (cm/s)/km");
fw.SetColumnName(6, "Rot, (cm/s)/km"); fw.SetColumnName(6, "Rot, (cm/s)/km");
fw.SetColumnName(7, "Okubo-Weiss parameter, (cm^2/s^2)/km^2"); 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); fw.SetParameters(pars);
for(size_t i = 0; i < data.N(); i++) 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.Div(i) == data.Fillval() ? NAN : data.Div(i));
fw.Write(data.Rot(i) == data.Fillval() ? NAN : data.Rot(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.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.Finalize();
fw.Close(); fw.Close();
@ -166,4 +113,4 @@ MString Data::ActionUV(const CLArgs& args)
} }
return ""; return "";
} };

81
include/basedata.h

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

16
include/data.h

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

13
include/layereddata.h

@ -1,25 +1,14 @@
#pragma once #pragma once
#include "DataAdapters/ncfilealt.h"
#include "ParseArgs.h"
#include "gsw.h" #include "gsw.h"
#include "mdatetime.h"
#include "ncfuncs.h" #include "ncfuncs.h"
#include "simple2ddata.h" #include "simple2ddata.h"
#include "uvdata.h"
#include <algorithm>
#include <memory> #include <memory>
#include <set>
using michlib::Ceil; using michlib::Ceil;
using michlib::DetGeoDomain; using michlib::DetGeoDomain;
using michlib::Floor; using michlib::Floor;
using michlib::GPL; using michlib::GPL;
using michlib::int2; using michlib::int2;
using michlib::int4;
using michlib::MDateTime;
using michlib::MString;
using michlib::NCFileA;
using michlib::ToGeoDomain;
class LayeredData: public NCFuncs class LayeredData: public NCFuncs
{ {
@ -148,8 +137,6 @@ class LayeredData: public NCFuncs
Data Read(const MString& vname, const BaseParameters* ip, size_t i) const; 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; } bool isOk() const { return nc.size() > 0; }
explicit operator bool() const { return nc.size() > 0; } explicit operator bool() const { return nc.size() > 0; }

1
include/ncfuncs.h

@ -1,6 +1,5 @@
#pragma once #pragma once
#include "DataAdapters/ncfilealt.h" #include "DataAdapters/ncfilealt.h"
#include "StringFunctions.h"
#include "mdatetime.h" #include "mdatetime.h"
#include <set> #include <set>
#include <tuple> #include <tuple>

186
include/odm.h

@ -1,180 +1,18 @@
#pragma once #pragma once
#include "AVISO.h" #include "actioninfo.h"
#include "AVISOLOCAL.h" #include "actiontsc.h"
#include "BFileW.h" #include "actionuv.h"
#include "HYCOM.h" #include "data.h"
#include "NEMO.h"
#include <utility>
#include <variant>
using michlib::BFileW;
using michlib::errmessage; using michlib::errmessage;
using michlib::GPL;
using michlib::message;
using DataVariants = std::variant<NEMOData, HYCOMData, AVISOData, AVISOLOCALData>; template<class T, ActionID id> consteval bool MustSupported()
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
{ {
public: static_assert(IsActionSupported<T, id>);
using TIndex = std::vector<size_t>; return true;
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 == "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
return "Unknown source: " + src;
return "";
}
MString ActionInfo(const CLArgs& args);
MString ActionTsc(const CLArgs& args);
MString ActionUV(const CLArgs& args);
}; };
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>();

74
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;

252
include/uvdata.h

@ -6,22 +6,35 @@ using michlib::M_PI;
using michlib::real; using michlib::real;
using michlib::Sqrt; 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 enum class Metric
{ {
EUCLID, EUCLID,
SPHERIC SPHERIC
}; };
class UVData class UVDataStorage
{ {
protected:
static constexpr real fillval = 1.0e10; 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; Metric metric;
UVDataStorage() = default;
real D(real lon1, real lat1, real lon2, real lat2) const real D(real lon1, real lat1, real lon2, real lat2) const
{ {
switch(metric) 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 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 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: public:
enum STPOINT enum STPOINT
{ {
@ -93,12 +69,6 @@ class UVData
STPOINT type = NOPOINT; 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& U(size_t i) const { return u[i]; }
const real& V(size_t i) const { return v[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); } 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& 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); } 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 N() const { return u.size(); }
size_t Nx() const { return nx; } size_t Nx() const { return nx; }
size_t Ny() const { return ny; } 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 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 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 Lon(size_t i) const { return Lon(i % nx, i / nx); }
real Lat(size_t i) const { return Lat(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; } void Resize(size_t sz)
bool IsCoast(size_t ix, size_t iy) const { return U(ix, iy) == fillval || V(ix, iy) == fillval; } {
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 Div(size_t i) const { return Div(i % nx, i / B::nx); }
real Rot(size_t i) const { return Rot(i % nx, i / nx); } real Rot(size_t i) const { return Rot(i % nx, i / B::nx); }
real Div(size_t ix, size_t iy) const real Div(size_t ix, size_t iy) const
{ {
if(!*this) return 0.0; if(!*this) return 0.0;
real ux = dUdX(ix, iy); real ux = dUdX(ix, iy);
real vy = dVdY(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 real Rot(size_t ix, size_t iy) const
{ {
if(!*this) return 0.0; if(!*this) return 0.0;
real vx = dVdX(ix, iy); real vx = dVdX(ix, iy);
real uy = dUdY(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 // Okubo-Weiss parameter
real OW(size_t i) const { return OW(i % nx, i / nx); } real OW(size_t i) const { return OW(i % nx, i / nx); }
@ -152,7 +272,7 @@ class UVData
real vx = dVdX(ix, iy); real vx = dVdX(ix, iy);
real vy = dVdY(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 sn = ux - vy;
real ss = vx + uy; real ss = vx + uy;
@ -162,7 +282,7 @@ class UVData
auto StablePoints(size_t ix, size_t iy) const 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(!*this) return points;
if(ix >= Nx() - 1 || iy >= Ny() - 1) 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; 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(des > 0.0)
{ {
if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0 && c > 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) else if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0)
return SADDLE; return UVDataStorage::SADDLE;
else else
return SKNOT; return UVDataStorage::SKNOT;
} }
else else
{ {
bool acyclcrit = (av * y + cv > 0.0); bool acyclcrit = (av * y + cv > 0.0);
if(b < 0.0) if(b < 0.0)
return acyclcrit ? SACICFOCUS : SCICFOCUS; return acyclcrit ? UVDataStorage::SACICFOCUS : UVDataStorage::SCICFOCUS;
else 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(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(x2 + ix, y2 + iy), LatR(x2 + ix, y2 + iy), PointType(x2, y2)); 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; return points;
} }
}; };

22
src/AVISOLOCAL.cpp

@ -234,25 +234,3 @@ template<class DataType> AVISOLOCALData::Data AVISOLOCALData::ReadVarRaw(const N
} }
return data; return data;
} }
UVData AVISOLOCALData::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), u.XStep(), u.YStep()};
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;
}

2
src/ParseArgs.cpp

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

19
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 "";
}

113
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 "";
}

1
src/basedata.cpp

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

39
src/data.cpp

@ -0,0 +1,39 @@
#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
return "Unknown source: " + src;
return "";
}

22
src/layereddata.cpp

@ -270,28 +270,6 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters*
return Data(); 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 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; real unitmul = 1.0;

4
src/ncfuncs.cpp

@ -95,8 +95,8 @@ std::tuple<MDateTime, time_t, bool> NCFuncs::Refdate(const MString& refdate)
if(*ci == "days") step = 3600 * 24; if(*ci == "days") step = 3600 * 24;
ci++; ci++;
} }
if(ci != words.end()) ci++; // skip "since" if(ci != words.end()) ci++; // skip "since"
if(ci != words.end()) rstr = *ci; // Day if(ci != words.end()) rstr = *ci; // Day
if(ci != words.end()) ci++; if(ci != words.end()) ci++;
if(ci != words.end()) rstr += " " + *ci; // Hours if(ci != words.end()) rstr += " " + *ci; // Hours
bool success = out.FromString(rstr); bool success = out.FromString(rstr);

39
src/odm.cpp

@ -4,7 +4,7 @@ inline void Usage(const MString& arg0)
{ {
message(arg0 + " (<key>=<value>...)"); message(arg0 + " (<key>=<value>...)");
message("Keys are:"); 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(" Keys for action=info. Print some information about dataset.");
message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL"); message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL");
message(" Keys for source=NEMO"); message(" Keys for source=NEMO");
@ -50,7 +50,7 @@ int main(int argc, char** argv)
} }
auto args = ParseArgs(argc, argv); auto args = ParseArgs(argc, argv);
MString action = args.contains("action") ? args["action"] : "tsc"; MString action = args.contains("action") ? args["action"] : "info";
Data data; Data data;
{ {
@ -61,19 +61,32 @@ int main(int argc, char** argv)
return 1; return 1;
} }
} }
Action act;
MString ret;
if(action == "info")
ret = data.ActionInfo(args);
else if(action == "tsc")
ret = data.ActionTsc(args);
else if(action == "uv")
ret = data.ActionUV(args);
else
{ {
errmessage("Unknown action " + action); auto ret = act.Init(args);
return 1; 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()) if(ret.Exist())
{ {
errmessage(ret); errmessage(ret);

Loading…
Cancel
Save