odm/actions/actionint.h

#pragma once
#include "BFileR.h"
#include "BFileW.h"
#include "actiondep.h"
#include "interpolation.h"
#include "ncfilew.h"
#include "ncfuncs.h"
#include <memory>

template<class T>
concept ReadIsInterpolable = requires {
 {
  std::declval<ReadType<T>>().GridPos(0.0, 0.0)
 } -> std::convertible_to<struct GridPoint>;
};

ADD_ACTION(INT, int, (ReadPSupported<Source> || ReadSupported<Source>)&&ReadIsInterpolable<Source>);

template<class D> MString ActionINT::DoAction(const CLArgs& args, D& ds)
{
 struct TPoint
 {
  struct Point2D p;
  MDateTime      t;
  real           lon0, lat0, dtime;
 };

 MDateTime refdate("1980-01-01");

 if(!args.contains("input")) return "No input file given";
 michlib::BFileR in;
 if(in.Open(args.at("input")) != ERR_NOERR) return "Can't open input file " + args.at("input");
 std::vector<TPoint> points;
 struct Region       reg;

 {
  michlib::uint loncol = in.Columns(), latcol = in.Columns(), timecol = in.Columns();
  michlib::uint lon0col = in.Columns(), lat0col = in.Columns(), dtimecol = in.Columns();
  for(uint i = 0; i < in.Columns(); i++)
  {
   if(in.ColumnName(i + 1) == "lon") loncol = i;
   if(in.ColumnName(i + 1) == "lat") latcol = i;
   if(in.ColumnName(i + 1) == "rtime") timecol = i;
   if(in.ColumnName(i + 1) == "lon0") lon0col = i;
   if(in.ColumnName(i + 1) == "lat0") lat0col = i;
   if(in.ColumnName(i + 1) == "time") dtimecol = i;
  }
  if(loncol >= in.Columns()) return "Lon column not found in input file " + args.at("input");
  if(latcol >= in.Columns()) return "Lat column not found in input file " + args.at("input");
  if(timecol >= in.Columns()) return "Time column not found in input file " + args.at("input");
  if(lon0col >= in.Columns()) return "Lon0 column not found in input file " + args.at("input");
  if(lat0col >= in.Columns()) return "Lat0 column not found in input file " + args.at("input");
  if(dtimecol >= in.Columns()) return "Age column not found in input file " + args.at("input");

  points.resize(in.Rows());
  reg.lonb = reg.lone = in[loncol][0];
  reg.latb = reg.late = in[latcol][0];
  for(uint i = 0; i < in.Rows(); i++)
  {
   points[i].p.x   = in[loncol][i];
   points[i].p.y   = in[latcol][i];
   points[i].t     = refdate + static_cast<time_t>(michlib::Round(in[timecol][i] * 86400));
   points[i].lon0  = in[lon0col][i];
   points[i].lat0  = in[lat0col][i];
   points[i].dtime = in[dtimecol][i];

   reg.lonb = std::min(reg.lonb, points[i].p.x);
   reg.lone = std::max(reg.lone, points[i].p.x);
   reg.latb = std::min(reg.latb, points[i].p.y);
   reg.late = std::max(reg.late, points[i].p.y);
  }
  std::ranges::sort(points, {}, &TPoint::t);
 }
 in.Close();

 auto resop = ds.Open(args);
 if(resop.Exist()) return "Can't open source: " + resop;

 michlib_internal::ParameterListEx pars;
 pars.UsePrefix("");
 pars.SetParameter("source", args.at("source"));
 pars.SetParameter("history", args.at("_cmdline"));

 MString varstring;
 if(args.contains("var"))
  varstring = args.at("var");
 else
 {
  if(args.contains("vars"))
   varstring = args.at("vars");
  else
  {
   if constexpr(HasDefVars<D>)
    varstring = ds.DefaultVars();
   else
    return "Variables not specified";
  }
 }

 auto                 vlist = michlib::Split_on_words(varstring, " ,", false);
 std::vector<MString> vnames{std::make_move_iterator(std::begin(vlist)), std::make_move_iterator(std::end(vlist))};
 {
  std::set<MString> used;
  for(const auto& vname: vnames)
  {
   if(used.contains(vname)) return "Duplicate variable " + vname + " in list " + varstring;
   if(ds.CheckVar(vname) == VarPresence::NONE) return "Variable " + vname + " not exists in this dataset";
   used.insert(vname);
  }
 }

 pars.SetParameter("variables", varstring);

 //int compress = 3;
 //if(args.contains("compress")) compress = args.at("compress").ToInt();

 std::unique_ptr<const BaseParameters> sourcepars;
 if constexpr(ParametersSupported<D>)
 {
  if constexpr(ParametersRequiredRegion<D>)
  {
   auto [p, err] = ds.Parameters(pars, args, reg);
   if(err.Exist()) return err;
   sourcepars.reset(p);
  }
  else
  {
   auto [p, err] = ds.Parameters(pars, args);
   if(err.Exist()) return err;
   sourcepars.reset(p);
  }
 }
 auto p = sourcepars.get();

 MString name = args.contains("out") ? args.at("out") : "out.bin";
 //MString outfmt = args.contains("format") ? args.at("format") : (GetExt(name) == "nc" ? "nc" : "bin");

 size_t          loind = ds.NTimes(), hiind = ds.NTimes(), curind = 0;
 michlib::BFileW fw;
 const MDateTime first = ds.Time(0), last = ds.Time(ds.NTimes() - 1);

 std::vector<LinearInterpolator<ReadType<D>>> lo, hi;

 {
  size_t ic = 0;
  fw.Create(name, 6 + vnames.size());
  fw.SetColumnName(++ic, "lon0");
  fw.SetColumnName(++ic, "lat0");
  fw.SetColumnName(++ic, "lon");
  fw.SetColumnName(++ic, "lat");
  fw.SetColumnName(++ic, "time");
  fw.SetColumnName(++ic, "rtime");
  for(size_t i = 0; i < vnames.size(); i++) fw.SetColumnName(++ic, vnames[i]);
  fw.SetParameters(pars);
 }

 for(size_t i = 0; i < points.size(); i++)
 {
  if(points[i].t < first || points[i].t > last) continue;
  while(!(points[i].t >= ds.Time(curind) && points[i].t <= ds.Time(curind + 1))) curind++;
  if(curind != loind && curind != hiind)
  {
   loind = curind;
   hiind = curind + 1;
   {
    auto temp = Read(ds, vnames, p, loind);
    if(temp.size() != vnames.size()) return "Can't read data";
    //lo.resize(temp.size());
    for(size_t j = 0; j < temp.size(); j++) lo.emplace_back(std::move(temp[j]));
   }
   {
    auto temp = Read(ds, vnames, p, hiind);
    if(temp.size() != vnames.size()) return "Can't read data";
    //hi.resize(temp.size());
    for(size_t j = 0; j < temp.size(); j++) hi.emplace_back(std::move(temp[j]));
   }
  }
  else if(curind == hiind)
  {
   loind = curind;
   hiind = curind + 1;
   lo    = std::move(hi);
   {
    auto temp = Read(ds, vnames, p, hiind);
    if(temp.size() != vnames.size()) return "Can't read data";
    //hi.resize(temp.size());
    for(size_t j = 0; j < temp.size(); j++) hi.emplace_back(std::move(temp[j]));
   }
  }
  auto step  = ds.Time(hiind) - ds.Time(loind);
  auto delta = points[i].t - ds.Time(loind);
  real trel  = static_cast<real>(delta) / static_cast<real>(step);

  fw.Write(points[i].lon0);
  fw.Write(points[i].lat0);
  fw.Write(points[i].p.x);
  fw.Write(points[i].p.y);
  fw.Write(points[i].dtime);
  fw.Write(DeltaDates(refdate, points[i].t));
  for(size_t iv = 0; iv < lo.size(); iv++)
  {
   real vlo = lo[iv](points[i].p);
   real vhi = hi[iv](points[i].p);
   fw.Write(vlo + (vhi - vlo) * trel);
  }
 }
 fw.Finalize();
 fw.Close();
 return "";
};