odm/actions/actionuv.h

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

using michlib::BFileW;

ADD_ACTION(UV, uv, ReadPSupported<Source> || ReadSupported<Source>);

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

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

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

 auto data = ReadUV(ds, p, tindexes);
 if(!data) return "Can't read data";

 // Main file
 MString name = args.contains("out") ? args.at("out") : "";
 if(name.Exist())
 {
  BFileW fw;
  fw.Create(name, 9);
  fw.SetColumnName(1, "Longitude");
  fw.SetColumnName(2, "Latitude");
  fw.SetColumnName(3, "u, cm/s");
  fw.SetColumnName(4, "v, cm/s");
  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++)
  {
   fw.Write(data.Lon(i));
   fw.Write(data.Lat(i));
   fw.Write(data.U(i) == data.Fillval() ? NAN : data.U(i));
   fw.Write(data.V(i) == data.Fillval() ? NAN : data.V(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.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();
 }

 // Filtered vectors file
 name = args.contains("velout") ? args.at("velout") : "";
 if(name.Exist())
 {
  size_t shiftx = args.contains("shiftx") ? args.at("shiftx").ToInteger<size_t>() : 0;
  size_t shifty = args.contains("shifty") ? args.at("shifty").ToInteger<size_t>() : 0;
  size_t skipx  = args.contains("skipx") ? args.at("skipx").ToInteger<size_t>() : 1;
  size_t skipy  = args.contains("skipy") ? args.at("skipy").ToInteger<size_t>() : 1;

  BFileW vel;
  vel.Create(name, 4);
  vel.SetColumnName(1, "Longitude");
  vel.SetColumnName(2, "Latitude");
  vel.SetColumnName(3, "u, cm/s");
  vel.SetColumnName(4, "v, cm/s");

  for(size_t ix = shiftx; ix < data.Nx(); ix += skipx)
   for(size_t iy = shifty; iy < data.Ny(); iy += skipy)
   {
    vel.Write(data.Lon(ix, iy));
    vel.Write(data.Lat(ix, iy));
    vel.Write(data.U(ix, iy) == data.Fillval() ? NAN : data.U(ix, iy));
    vel.Write(data.V(ix, iy) == data.Fillval() ? NAN : data.V(ix, iy));
   }
  vel.Finalize();
  vel.Close();
 }

 // Stationary points
 name = args.contains("stpout") ? args.at("stpout") : "";
 if(name.Exist())
 {
  BFileW stp;
  stp.Create(name, 3);
  stp.SetColumnName(1, "Longitude");
  stp.SetColumnName(2, "Latitude");
  stp.SetColumnName(3, "Stability (0 - saddle, 1 - st. anticicl. focus, 2 - st. knot, 3 - unst. anticicl. focus, 4 - unst. knot, 5 - st. cicl. focus, 6 - unst. cicl. focus)");

  for(size_t ix = 0; ix < data.Nx() - 1; ix++)
   for(size_t iy = 0; iy < data.Ny() - 1; iy++)
   {
    auto points = data.StablePoints(ix, iy);
    for(size_t i = 0; i < points.size(); i++)
    {
     stp.Write(points[i].x);
     stp.Write(points[i].y);
     stp.Write(michlib::int_cast<size_t>(points[i].type));
    }
   }
  stp.Finalize();
  stp.Close();
 }

 return "";
};