odm/actions/actiongenintfile.h

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

using michlib::BFileW;
using michlib::Cos;
using michlib::M_PI;

ADD_ACTION(GenIntFile, genintfile, ReadIsGrid<Source>);

template<class D> MString ActionGenIntFile::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"));

 if(!args.contains("out")) return "Output file name not specified";

 VelSource mode = VelSource::AUTOSEL;
 if(args.contains("geostrophic"))
 {
  auto modestr = args.at("geostrophic");
  if(modestr == "" || modestr == "ssh")
   mode = VelSource::GEOSTROPHICSSH;
  else if(modestr == "surf" || modestr == "surface")
   mode = VelSource::GEOSTROPHIC;
  else if(modestr == "nongeo")
   mode = VelSource::STANDART;
  else
   return "Parameter \"geostrophic\" have incorrect value " + modestr;
 }

 auto [tindexes, err] = GetTIndexes(ds, args, pars);
 if(err.Exist()) return err;
 if(tindexes.size() < 10) return "Too few time moments";
 // Check isotropy
 for(size_t it = 0; it < tindexes.size() - 1; it++)
  if(ds.Time(tindexes[it + 1]) - ds.Time(tindexes[it]) != ds.Time(tindexes[1]) - ds.Time(tindexes[0])) return "The time step is not isotropic";

 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.at("out");
 BFileW  fw;
 fw.Create(name, 2);
 fw.SetColumnName(1, "u");
 fw.SetColumnName(2, "v");

 for(size_t it = 0; it < tindexes.size(); it++)
 {
  auto data = ReadUV(ds, p, tindexes[it], mode);
  if(!data) return "Can't read data";
  if(it == 0)
  {
   real lonb = data.Lon(0, 0), latb = data.Lat(0, 0), lone = data.Lon(data.Nx() - 1, data.Ny() - 1), late = data.Lat(data.Nx() - 1, data.Ny() - 1);
   fw.UsePrefix("");
   fw.SetParameter("nx", data.Nx());
   fw.SetParameter("ny", data.Ny());
   fw.SetParameter("dx", data.XStep() * 60.0);
   fw.SetParameter("dy", data.YStep() * 60.0);
   fw.SetParameter("nt", tindexes.size());
   fw.SetParameter("dt", (ds.Time(tindexes[1]) - ds.Time(tindexes[0])).D());
   fw.SetParameter("FillValue", data.Fillval());
   fw.UsePrefix("Info");
   switch(mode)
   {
   case(VelSource::AUTOSEL):
   {
    fw.SetParameter("Mode", "autoselect between standart and geostrophic");
    break;
   }
   case(VelSource::STANDART):
   {
    fw.SetParameter("Mode", "standart");
    break;
   }
   case(VelSource::GEOSTROPHIC):
   {
    fw.SetParameter("Mode", "geostrophic");
    break;
   }
   case(VelSource::GEOSTROPHICSSH):
   {
    fw.SetParameter("Mode", "geostrophic from ssh");
    break;
   }
   }
   fw.SetParameter("lonb", lonb);
   fw.SetParameter("latb", latb);
   fw.SetParameter("lone", lone);
   fw.SetParameter("late", late);
   fw.SetParameter("DataSource", ds.DataTitle());
   fw.SetParameter("xy2lon", MString(lonb) + "+%x/60");
   fw.SetParameter("xy2lon_rp", MString("POP 60 DIV ") + lonb + " ADD");
   fw.SetParameter("lonlat2x", MString("(%lon-") + lonb + ")*60");
   fw.SetParameter("xy2lat", MString(latb) + "+%y/60");
   fw.SetParameter("xy2lat_rp", MString("EXCH POP 60 DIV ") + latb + " ADD");
   fw.SetParameter("lonlat2y", MString("(%lat-") + latb + ")*60");
   fw.SetParameter("lonlatuv2dotx", "%u*(0.864/1.852)/Cos(%lat*M_PI/180)");
   fw.SetParameter("lonlatuv2doty", "%v*(0.864/1.852)");
   fw.SetParameter("dotxdotylonlat2u_rp", "EXCH POP EXCH POP COSD 1.852 0.864 DIV MUL MUL");
   fw.SetParameter("dotxdotylonlat2v_rp", "POP POP 1.852 0.864 DIV MUL EXCH POP");
   fw.SetParameter("dotxdotylonlat2u", "(1.852/0.864)*Cos(%lat*M_PI/180)*%dotx");
   fw.SetParameter("dotxdotylonlat2v", "(1.852/0.864)*%doty");
   fw.SetParameter("BeginDate", ds.Time(tindexes.front()).ToTString());
   fw.SetParameter("EndDate", ds.Time(tindexes.back()).ToTString());
   fw.SetParameter("Timestep", (ds.Time(tindexes[1]) - ds.Time(tindexes[0])).Seconds());
   fw.SetParameter("DataID", michlib::UniqueId());
   fw.SetParameter("Creation command", args.at("_cmdline"));
  }

  for(size_t ilat = 0; ilat < data.Ny(); ilat++)
  {
   for(size_t ilon = 0; ilon < data.Nx(); ilon++)
   {
    if(data.IsCoast(ilon, ilat))
    {
     fw.Write(data.Fillval());
     fw.Write(data.Fillval());
    }
    else
    {
     fw.Write(data.U(ilon, ilat) * (0.864 / 1.852) / Cos(data.Lat(ilon, ilat) * M_PI / 180.0));
     fw.Write(data.V(ilon, ilat) * (0.864 / 1.852));
    }
   }
  }
 }
 fw.Finalize();
 fw.Close();

 return "";
};