odm/actions/actiontsc.h

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

using michlib::BFileW;

ADD_ACTION(TSC, tsc, ReadPSupported<Source> || ReadSupported<Source>);

template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
{
 auto [reg, regerr] = GetRegion<D>(args);
 if(regerr.Exist()) return regerr;

 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"));

 auto [tindexes, err] = GetTIndexes(ds, args, pars);
 if(err.Exist()) return err;

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

 bool average  = args.contains("average");
 bool gradient = args.contains("gradient");

 if(gradient)
  if constexpr(!ReadIs2DGeoArray<D>) return "Gradient calculation not supported for this source";

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

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

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

 TimeData tdata(ds, tindexes);

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

 if(outfmt == "bin" && !average && tindexes.size() > 1) return "Multiple time moments does'nt supported by the bin format";

 bool    headwrited = false;
 NCFileW ncfw;

 for(size_t it = 0; it < tindexes.size(); it++)
 {
  auto data = average ? Read(ds, vnames, p, tindexes) : Read(ds, vnames, p, tindexes[it]);
  if(data.size() != vnames.size()) return "Can't read data";
  if(!headwrited)
   for(size_t i = 0; i < vnames.size(); i++)
   {
    if(!data[i].Unit().Exist()) michlib::errmessage("Unknown measurement unit for variable " + vnames[i] + "!");
    if(!data[i].StandartName().Exist()) data[i].SetStandartName(NCFuncs::Name2StName(vnames[i]));
    if(!data[i].LongName().Exist()) data[i].SetLongName(NCFuncs::Name2LongName(vnames[i]));
   }

  if(outfmt == "bin")
  {
   BFileW fw;

   size_t ic = 0;
   fw.Create(name, 2 + data.size());
   fw.SetColumnName(++ic, "Longitude");
   fw.SetColumnName(++ic, "Latitude");
   for(size_t i = 0; i < data.size(); i++)
    fw.SetColumnName(++ic, (gradient ? "gradient of " : "") + vnames[i] + ", " + (data[i].Unit().Exist() ? data[i].Unit() : "unknown") + (gradient ? "/km" : ""));
   fw.SetParameters(pars);
   for(size_t r = 0; r < data[0].N(); r++)
   {
    fw.Write(data[0].Lon(r));
    fw.Write(data[0].Lat(r));
    if(gradient)
    {
     if constexpr(ReadIs2DGeoArray<D>)
      for(size_t i = 0; i < data.size(); i++)
      {
       auto val = data[i].Grad(r);
       fw.Write(val == data[i].Fillval() ? NAN : val);
      }
    }
    else
     for(size_t i = 0; i < data.size(); i++) fw.Write(data[i].IsFill(r) ? NAN : data[i](r));
   }
   fw.Finalize();
   fw.Close();
  }
  else if(outfmt == "nc" || outfmt == "netcdf")
  {
   MString err;

   if(!err.Exist() && !headwrited) err = ncfw.Create(data[0], name, compress);
   if(!err.Exist() && !headwrited) err = ncfw.AddTimeData(tdata, !average);
   if(!err.Exist() && !headwrited) err = ncfw.AddAtts(pars);
   for(size_t i = 0; i < data.size(); i++)
    if(!err.Exist() && !headwrited)
     err = ncfw.AddVariable((gradient ? "G" : "") + vnames[i], data[i].StandartName(), data[i].LongName(), data[i].Unit() + (gradient ? "/km" : ""), data[i].Comment());
   if(!err.Exist() && !headwrited) err = ncfw.WriteGrid(data[0]);
   for(size_t i = 0; i < data.size(); i++)
    if(!err.Exist())
    {
     if(gradient)
     {
      if constexpr(ReadIs2DGeoArray<D>)
       err = average ? ncfw.WriteVariable(data[i], "G" + vnames[i], [&data = std::as_const(data[i])](size_t i, size_t j) { return data.Grad(i, j); })
                     : ncfw.WriteVariable(data[i], "G" + vnames[i], [&data = std::as_const(data[i])](size_t i, size_t j) { return data.Grad(i, j); }, it);
     }
     else
      err = average ? ncfw.WriteVariable(data[i], vnames[i]) : ncfw.WriteVariable(data[i], vnames[i], it);
    }
   if(err.Exist()) return err;
  }
  else
   return "Unknown format: " + outfmt;

  if(average) break;
  headwrited = true;
 }
 ncfw.Close();

 return "";
};