Support for time column in netcdf files generated by uv and tsc actions

actions/actiontsc.h

@@ -41,6 +41,11 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
   }
  }
 
+ bool average = args.contains("average");
+
+ 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))};
  {
@@ -73,60 +78,68 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
  }
  auto p = sourcepars.get();
 
- auto data = Read(ds, vnames, p, tindexes);
- if(data.size() != vnames.size()) return "Can't read data";
-
- 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]));
- }
+ 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")
+ 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++)
  {
-  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, vnames[i] + ", " + (data[i].Unit().Exist() ? data[i].Unit() : "unknown"));
-  fw.SetParameters(pars);
-  for(size_t r = 0; r < data[0].N(); r++)
+  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")
   {
-   fw.Write(data[0].Lon(r));
-   fw.Write(data[0].Lat(r));
-   for(size_t i = 0; i < data.size(); i++) fw.Write(data[i].IsFill(r) ? NAN : data[i](r));
+   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, vnames[i] + ", " + (data[i].Unit().Exist() ? data[i].Unit() : "unknown"));
+   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));
+    for(size_t i = 0; i < data.size(); i++) fw.Write(data[i].IsFill(r) ? NAN : data[i](r));
+   }
+   fw.Finalize();
+   fw.Close();
   }
-  fw.Finalize();
-  fw.Close();
-  return "";
- }
+  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(vnames[i], data[i].StandartName(), data[i].LongName(), data[i].Unit(), data[i].Comment());
+   if(!err.Exist() && !headwrited) err = ncfw.WriteGrid(data[0]);
+   for(size_t i = 0; i < data.size(); i++)
+    if(!err.Exist()) 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(outfmt == "nc" || outfmt == "netcdf")
- {
-  int     compress = 3;
-  NCFileW fw;
-  MString err;
-
-  if(args.contains("compress")) compress = args.at("compress").ToInt();
-
-  if(!err.Exist()) err = fw.Create(data[0], name, compress);
-  if(!err.Exist()) err = fw.AddAtts(pars);
-  for(size_t i = 0; i < data.size(); i++)
-   if(!err.Exist()) err = fw.AddVariable(vnames[i], data[i].StandartName(), data[i].LongName(), data[i].Unit(), data[i].Comment());
-  if(!err.Exist()) err = fw.WriteGrid(data[0]);
-  for(size_t i = 0; i < data.size(); i++)
-   if(!err.Exist()) err = fw.WriteVariable(data[i], vnames[i]);
-  if(err.Exist()) return err;
-  fw.Close();
-
-  return "";
+  if(average) break;
+  headwrited = true;
  }
+ ncfw.Close();
 
- return "Unknown format: " + outfmt;
+ return "";
 };

actions/actionuv.cpp

@@ -20,9 +20,9 @@ void UVMethods::StPoints::WriteBinBile(const MString& name, const michlib_intern
  stp.Close();
 }
 
-MString UVMethods::StPoints::WriteNcFile(const MString& name, const michlib_internal::ParameterListEx& pars, const MString& history, int comp) const
+MString UVMethods::StPoints::CreateNcFile(const MString& name, const michlib_internal::ParameterListEx& pars, const MString& history, int comp, const TimeData& tdata, bool timedep)
 {
- NCFileWBase nc;
+ tdep = timedep;
 
  const MString xname = lonlat ? "longitude" : "x";
  const MString yname = lonlat ? "latitude" : "y";
@@ -33,8 +33,24 @@ MString UVMethods::StPoints::WriteNcFile(const MString& name, const michlib_inte
  nc.AddAtt("history", history);
  nc.AddAtts(pars);
  nc.AddDim("i", N());
- nc.AddVar(xname, NC_FLOAT, "i");
- nc.AddVar(yname, NC_FLOAT, "i");
+ nc.AddDim("time", tdata.steps.size());
+
+ nc.AddVar("time", decltype(nc)::Type2NCType<decltype(tdata.steps)::value_type>, "time");
+ nc.SetComp("time", comp);
+ nc.AddAtt("time", "standard_name", "time");
+ nc.AddAtt("time", "long_name", "time");
+ nc.AddAtt("time", "units", tdata.UnitName());
+
+ if(tdep)
+ {
+  nc.AddVar(xname, NC_FLOAT, "time", "i");
+  nc.AddVar(yname, NC_FLOAT, "time", "i");
+ }
+ else
+ {
+  nc.AddVar(xname, NC_FLOAT, "i");
+  nc.AddVar(yname, NC_FLOAT, "i");
+ }
  nc.SetComp(xname, comp);
  nc.SetComp(yname, comp);
 
@@ -51,20 +67,38 @@ MString UVMethods::StPoints::WriteNcFile(const MString& name, const michlib_inte
   nc.AddAtt(yname, "long_name", "y-coordinate");
  }
 
- nc.AddVar("type", NC_UBYTE, "i");
+ if(tdep)
+  nc.AddVar("type", NC_UBYTE, "time", "i");
+ else
+  nc.AddVar("type", NC_UBYTE, "i");
  nc.SetComp("type", comp);
  nc.AddAtt("type", "long_name",
            "Stationary point type, 0 - saddle, 1 - st. anticicl. focus, 2 - st. knot, 3 - unst. anticicl. focus, 4 - unst. knot, 5 - st. cicl. focus, 6 - unst. cicl. focus");
 
- if(!nc) return "Can't set grid in the netcdf file " + name + ": " + nc.ErrMessage();
+ nc.WriteVar("time", tdata.steps.data());
 
- nc.EndDef();
+ if(!nc) return "Can't set grid in the netcdf file " + name + ": " + nc.ErrMessage();
+ return "";
+}
 
- nc.WriteVar(xname, x.data());
- nc.WriteVar(yname, y.data());
- nc.WriteVar("type", t.data());
+MString UVMethods::StPoints::WriteNcFile(size_t it)
+{
+ const MString xname = lonlat ? "longitude" : "x";
+ const MString yname = lonlat ? "latitude" : "y";
 
- if(!nc) return "Can't write data to the netcdf file " + name + ": " + nc.ErrMessage();
+ if(tdep)
+ {
+  nc.WriteVar(xname, it, x.data());
+  nc.WriteVar(yname, it, y.data());
+  nc.WriteVar("type", it, t.data());
+ }
+ else
+ {
+  nc.WriteVar(xname, x.data());
+  nc.WriteVar(yname, y.data());
+  nc.WriteVar("type", t.data());
+ }
 
+ if(!nc) return MString("Can't write data to the netcdf file: ") + nc.ErrMessage();
  return "";
 }

actions/actionuv.h

@@ -70,6 +70,8 @@ class UVMethods
   std::vector<real>           x, y;
   std::vector<michlib::uint1> t;
   bool                        lonlat;
+  NCFileWBase                 nc;
+  bool                        tdep;
 
   size_t N() const { return t.size(); }
 
@@ -86,9 +88,19 @@ class UVMethods
    }
   }
 
+  void Reset()
+  {
+   x.clear();
+   y.clear();
+   t.clear();
+  }
+
   void WriteBinBile(const MString& name, const michlib_internal::ParameterListEx& pars) const;
 
-  MString WriteNcFile(const MString& name, const michlib_internal::ParameterListEx& pars, const MString& history, int comp) const;
+  MString CreateNcFile(const MString& name, const michlib_internal::ParameterListEx& pars, const MString& history, int comp, const TimeData& tdata, bool timedep);
+  MString WriteNcFile(size_t it);
+
+  void CloseNcFile() { nc.Close(); }
  };
 };
 
@@ -112,6 +124,11 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
  auto [tindexes, err] = GetTIndexes(ds, args, pars);
  if(err.Exist()) return err;
 
+ bool average = args.contains("average");
+
+ int compress = 3;
+ if(args.contains("compress")) compress = args.at("compress").ToInt();
+
  std::unique_ptr<const BaseParameters> sourcepars;
  if constexpr(ParametersSupported<D>)
  {
@@ -130,160 +147,188 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
  }
  auto p = sourcepars.get();
 
- auto data = ReadUV(ds, p, tindexes, mode);
- if(!data) return "Can't read data";
+ TimeData tdata(ds, tindexes);
 
- // Main file
- MString name   = args.contains("out") ? args.at("out") : "";
- MString outfmt = args.contains("outformat") ? args.at("outformat") : (GetExt(name) == "nc" ? "nc" : "bin");
+ MString name      = args.contains("out") ? args.at("out") : "";
+ MString outfmt    = args.contains("outformat") ? args.at("outformat") : (GetExt(name) == "nc" ? "nc" : "bin");
+ MString namevel   = args.contains("velout") ? args.at("velout") : "";
+ MString outfmtvel = args.contains("veloutformat") ? args.at("veloutformat") : (GetExt(name) == "nc" ? "nc" : "bin");
+ MString namestp   = args.contains("stpout") ? args.at("stpout") : "";
+ MString outfmtstp = args.contains("stpoutformat") ? args.at("stpoutformat") : (GetExt(name) == "nc" ? "nc" : "bin");
 
- MString u = data.Unit().Exist() ? data.Unit() : "unknown", d = data.DUnit().Exist() ? data.DUnit() : "unknown";
+ 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;
 
- if(name.Exist())
- {
-  if(outfmt == "bin")
-  {
-   BFileW fw;
-   fw.Create(name, 9);
-   fw.SetColumnName(1, "Longitude");
-   fw.SetColumnName(2, "Latitude");
-   fw.SetColumnName(3, "u, " + u);
-   fw.SetColumnName(4, "v, " + u);
-   fw.SetColumnName(5, "Div, (" + u + ")/" + d);
-   fw.SetColumnName(6, "Rot, (" + u + ")/" + d);
-   fw.SetColumnName(7, "Okubo-Weiss parameter, (" + u + ")2/" + d + "2");
-   fw.SetColumnName(8, "Kinetic energy, (" + u + ")2");
-   fw.SetColumnName(9, "Eddy kinetic energy, (" + u + ")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();
-  }
-  else if(outfmt == "nc" || outfmt == "netcdf")
-  {
-   int     compress = 3;
-   NCFileW fw;
-   MString err;
-
-   if(args.contains("compress")) compress = args.at("compress").ToInt();
-
-   if(!err.Exist()) err = fw.Create(data, name, compress);
-   if(!err.Exist()) err = fw.AddAtts(pars);
-   if(!err.Exist()) err = fw.AddVariable("u", "", "Eastward velocity", u, "");
-   if(!err.Exist()) err = fw.AddVariable("v", "", "Northward velocity", u, "");
-   if(!err.Exist()) err = fw.AddVariable("div", "", "Velocity divergence", "(" + u + ")/" + d, "");
-   if(!err.Exist()) err = fw.AddVariable("rot", "", "Velocity rotor", "(" + u + ")/" + d, "");
-   if(!err.Exist()) err = fw.AddVariable("ow", "", "Okubo-Weiss parameter", "(" + u + ")2/" + d + "2", "");
-   if(!err.Exist()) err = fw.AddVariable("ke", "", "Squared velocity module, u^2+v^2", "(" + u + ")2", "");
-   if(!err.Exist()) err = fw.AddVariable("eke", "", "Squared velocity dispersion aka eddy kinetic energy, <u^2+v^2>-<u>^2-<v>^2", "(" + u + ")2", "");
-
-   if(!err.Exist()) err = fw.WriteGrid(data);
-
-   if(!err.Exist()) err = fw.WriteVariable(data, "u", [&data = std::as_const(data)](size_t i, size_t j) { return data.U(i, j); });
-   if(!err.Exist()) err = fw.WriteVariable(data, "v", [&data = std::as_const(data)](size_t i, size_t j) { return data.V(i, j); });
-   if(!err.Exist()) err = fw.WriteVariable(data, "div", [&data = std::as_const(data)](size_t i, size_t j) { return data.Div(i, j); });
-   if(!err.Exist()) err = fw.WriteVariable(data, "rot", [&data = std::as_const(data)](size_t i, size_t j) { return data.Rot(i, j); });
-   if(!err.Exist()) err = fw.WriteVariable(data, "ow", [&data = std::as_const(data)](size_t i, size_t j) { return data.OW(i, j); });
-   if(!err.Exist()) err = fw.WriteVariable(data, "ke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j); });
-   if(!err.Exist())
-    err = fw.WriteVariable(data, "eke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j) - (data.U(i, j) * data.U(i, j) + data.V(i, j) * data.V(i, j)); });
-   if(err.Exist()) return err;
+ if(!average && tindexes.size() > 1 && ((name.Exist() && outfmt == "bin") || (namevel.Exist() && outfmtvel == "bin") || (namestp.Exist() && outfmtstp == "bin")))
+  return "Multiple time moments does'nt supported by the bin format";
 
-   fw.Close();
-  }
-  else
-   return "Unknown format: " + outfmt;
- }
+ bool headwrited = false;
+
+ MString velunit, distunit;
+ NCFileW fw, fwfilt;
+
+ StPoints stp(true);
 
- // Filtered vectors file
- name   = args.contains("velout") ? args.at("velout") : "";
- outfmt = args.contains("veloutformat") ? args.at("veloutformat") : (GetExt(name) == "nc" ? "nc" : "bin");
- if(name.Exist())
+ for(size_t it = 0; it < tindexes.size(); it++)
  {
-  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;
+  auto data = average ? ReadUV(ds, p, tindexes, mode) : ReadUV(ds, p, tindexes[it], mode);
+  if(!headwrited)
+  {
+   velunit  = data.Unit().Exist() ? data.Unit() : "unknown";
+   distunit = data.DUnit().Exist() ? data.DUnit() : "unknown";
+  }
 
-  Sparser<decltype(data)> sdata(data, shiftx, shifty, skipx, skipy);
-  if(outfmt == "bin")
+  // Main file
+  if(name.Exist())
   {
-   BFileW vel;
-   vel.Create(name, 4);
-   vel.SetColumnName(1, "Longitude");
-   vel.SetColumnName(2, "Latitude");
-   vel.SetColumnName(3, "u, " + u);
-   vel.SetColumnName(4, "v, " + u);
-   vel.SetParameters(pars);
-
-   for(size_t ix = 0; ix < sdata.Nx(); ix++)
-    for(size_t iy = 0; iy < sdata.Ny(); iy++)
+   if(outfmt == "bin")
+   {
+    BFileW fw;
+    fw.Create(name, 9);
+    fw.SetColumnName(1, "Longitude");
+    fw.SetColumnName(2, "Latitude");
+    fw.SetColumnName(3, "u, " + velunit);
+    fw.SetColumnName(4, "v, " + velunit);
+    fw.SetColumnName(5, "Div, (" + velunit + ")/" + distunit);
+    fw.SetColumnName(6, "Rot, (" + velunit + ")/" + distunit);
+    fw.SetColumnName(7, "Okubo-Weiss parameter, (" + velunit + ")2/" + distunit + "2");
+    fw.SetColumnName(8, "Kinetic energy, (" + velunit + ")2");
+    fw.SetColumnName(9, "Eddy kinetic energy, (" + velunit + ")2");
+    fw.SetParameters(pars);
+
+    for(size_t i = 0; i < data.N(); i++)
     {
-     vel.Write(sdata.Lon(ix, iy));
-     vel.Write(sdata.Lat(ix, iy));
-     vel.Write(sdata.U(ix, iy) == sdata.Fillval() ? NAN : sdata.U(ix, iy));
-     vel.Write(sdata.V(ix, iy) == sdata.Fillval() ? NAN : sdata.V(ix, iy));
+     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))));
     }
-   vel.Finalize();
-   vel.Close();
+    fw.Finalize();
+    fw.Close();
+   }
+   else if(outfmt == "nc" || outfmt == "netcdf")
+   {
+    MString err;
+
+    if(!err.Exist() && !headwrited) err = fw.Create(data, name, compress);
+    if(!err.Exist() && !headwrited) err = fw.AddTimeData(tdata, !average);
+    if(!err.Exist() && !headwrited) err = fw.AddAtts(pars);
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("u", "", "Eastward velocity", velunit, "");
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("v", "", "Northward velocity", velunit, "");
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("div", "", "Velocity divergence", "(" + velunit + ")/" + distunit, "");
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("rot", "", "Velocity rotor", "(" + velunit + ")/" + distunit, "");
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("ow", "", "Okubo-Weiss parameter", "(" + velunit + ")2/" + distunit + "2", "");
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("ke", "", "Squared velocity module, u^2+v^2", "(" + velunit + ")2", "");
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("eke", "", "Squared velocity dispersion aka eddy kinetic energy, <u^2+v^2>-<u>^2-<v>^2", "(" + velunit + ")2", "");
+    if(!err.Exist() && !headwrited) err = fw.WriteGrid(data);
+
+    if(!err.Exist())
+     err = fw.WriteVariable(
+         data, "u", [&data = std::as_const(data)](size_t i, size_t j) { return data.U(i, j); }, it);
+    if(!err.Exist())
+     err = fw.WriteVariable(
+         data, "v", [&data = std::as_const(data)](size_t i, size_t j) { return data.V(i, j); }, it);
+    if(!err.Exist())
+     err = fw.WriteVariable(
+         data, "div", [&data = std::as_const(data)](size_t i, size_t j) { return data.Div(i, j); }, it);
+    if(!err.Exist())
+     err = fw.WriteVariable(
+         data, "rot", [&data = std::as_const(data)](size_t i, size_t j) { return data.Rot(i, j); }, it);
+    if(!err.Exist())
+     err = fw.WriteVariable(
+         data, "ow", [&data = std::as_const(data)](size_t i, size_t j) { return data.OW(i, j); }, it);
+    if(!err.Exist())
+     err = fw.WriteVariable(
+         data, "ke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j); }, it);
+    if(!err.Exist())
+     err = fw.WriteVariable(
+         data, "eke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j) - (data.U(i, j) * data.U(i, j) + data.V(i, j) * data.V(i, j)); }, it);
+    if(err.Exist()) return err;
+   }
+   else
+    return "Unknown format: " + outfmt;
   }
-  else if(outfmt == "nc" || outfmt == "netcdf")
-  {
-   int     compress = 3;
-   NCFileW fw;
-   MString err;
-
-   if(args.contains("compress")) compress = args.at("compress").ToInt();
-
-   if(!err.Exist()) err = fw.Create(sdata, name, compress);
-   if(!err.Exist()) err = fw.AddAtts(pars);
-   if(!err.Exist()) err = fw.AddVariable("u", "", "Eastward velocity", u, "");
-   if(!err.Exist()) err = fw.AddVariable("v", "", "Northward velocity", u, "");
-
-   if(!err.Exist()) err = fw.WriteGrid(sdata);
-
-   if(!err.Exist()) err = fw.WriteVariable(sdata, "u", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.U(i, j); });
-   if(!err.Exist()) err = fw.WriteVariable(sdata, "v", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.V(i, j); });
-   if(err.Exist()) return err;
 
-   fw.Close();
+  // Filtered file
+  if(namevel.Exist())
+  {
+   Sparser<decltype(data)> sdata(data, shiftx, shifty, skipx, skipy);
+   if(outfmtvel == "bin")
+   {
+    BFileW vel;
+    vel.Create(name, 4);
+    vel.SetColumnName(1, "Longitude");
+    vel.SetColumnName(2, "Latitude");
+    vel.SetColumnName(3, "u, " + velunit);
+    vel.SetColumnName(4, "v, " + velunit);
+    vel.SetParameters(pars);
+
+    for(size_t ix = 0; ix < sdata.Nx(); ix++)
+     for(size_t iy = 0; iy < sdata.Ny(); iy++)
+     {
+      vel.Write(sdata.Lon(ix, iy));
+      vel.Write(sdata.Lat(ix, iy));
+      vel.Write(sdata.U(ix, iy) == sdata.Fillval() ? NAN : sdata.U(ix, iy));
+      vel.Write(sdata.V(ix, iy) == sdata.Fillval() ? NAN : sdata.V(ix, iy));
+     }
+    vel.Finalize();
+    vel.Close();
+   }
+   else if(outfmtvel == "nc" || outfmtvel == "netcdf")
+   {
+    MString err;
+
+    if(!err.Exist() && !headwrited) err = fwfilt.Create(sdata, name, compress);
+    if(!err.Exist() && !headwrited) err = fwfilt.AddTimeData(tdata, !average);
+    if(!err.Exist() && !headwrited) err = fwfilt.AddAtts(pars);
+    if(!err.Exist() && !headwrited) err = fwfilt.AddVariable("u", "", "Eastward velocity", velunit, "");
+    if(!err.Exist() && !headwrited) err = fwfilt.AddVariable("v", "", "Northward velocity", velunit, "");
+    if(!err.Exist() && !headwrited) err = fwfilt.WriteGrid(sdata);
+
+    if(!err.Exist())
+     err = fwfilt.WriteVariable(
+         sdata, "u", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.U(i, j); }, it);
+    if(!err.Exist())
+     err = fwfilt.WriteVariable(
+         sdata, "v", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.V(i, j); }, it);
+    if(err.Exist()) return err;
+   }
+   else
+    return "Unknown format: " + outfmtvel;
   }
-  else
-   return "Unknown format: " + outfmt;
- }
-
- // Stationary points
- name   = args.contains("stpout") ? args.at("stpout") : "";
- outfmt = args.contains("stpoutformat") ? args.at("stpoutformat") : (GetExt(name) == "nc" ? "nc" : "bin");
- if(name.Exist())
- {
-  StPoints p(true);
-  for(size_t ix = 0; ix < data.Nx() - 1; ix++)
-   for(size_t iy = 0; iy < data.Ny() - 1; iy++) p.Add(data.StablePoints(ix, iy));
 
-  if(outfmt == "bin")
-   p.WriteBinBile(name, pars);
-  else if(outfmt == "nc" || outfmt == "netcdf")
+  // Stationary points
+  if(namestp.Exist())
   {
-   int     compress = args.contains("compress") ? args.at("compress").ToInt() : 3;
-   MString err      = p.WriteNcFile(name, pars, args.at("_cmdline"), compress);
-   if(err.Exist()) return err;
+   stp.Reset();
+   for(size_t ix = 0; ix < data.Nx() - 1; ix++)
+    for(size_t iy = 0; iy < data.Ny() - 1; iy++) stp.Add(data.StablePoints(ix, iy));
+
+   if(outfmtstp == "bin")
+    stp.WriteBinBile(name, pars);
+   else if(outfmtstp == "nc" || outfmtstp == "netcdf")
+   {
+    MString err;
+    if(!err.Exist() && !headwrited) err = stp.CreateNcFile(name, pars, args.at("_cmdline"), compress, tdata, !average);
+    if(!err.Exist()) err = stp.WriteNcFile(it);
+    if(err.Exist()) return err;
+   }
+   else
+    return "Unknown format: " + outfmt;
   }
-  else
-   return "Unknown format: " + outfmt;
- }
 
+  if(average) break;
+  headwrited = true;
+ }
+ fw.Close();
+ fwfilt.Close();
+ stp.CloseNcFile();
  return "";
 };

include/ncfilew.h

@@ -1,5 +1,6 @@
 #pragma once
 #include "MString.h"
+#include "actiondep.h"
 #include "traits.h"
 #include <algorithm>
 #include <netcdf.h>
@@ -366,6 +367,7 @@ class NCFileW: public NCFileWBase
 
  Type type = UNKNOWN;
  int  compress;
+ bool tdep = false;
 
  template<class D> static constexpr Type DetType()
  {
@@ -403,9 +405,12 @@ class NCFileW: public NCFileWBase
   type = UNKNOWN;
  }
 
+ MString AddTimeData(const TimeData& tdata, bool tisindex);
  MString AddVariable(const MString& name, const MString& stname, const MString& lname, const MString& units, const MString& comment);
 
- template<class D, class Op> MString WriteVariable(const D& data, const MString& name, Op op)
+ template<class D, class Op>
+ requires(std::is_invocable_v<Op, size_t> || std::is_invocable_v<Op, size_t, size_t>)
+ MString WriteVariable(const D& data, const MString& name, Op op, size_t tind)
  {
   static constexpr auto dtype = DetType<D>();
   if(type == UNKNOWN) return "File not open";
@@ -419,7 +424,10 @@ class NCFileW: public NCFileWBase
    const size_t c = data.N();
    float        buf[c];
    for(size_t ix = 0; ix < c; ix++) buf[ix] = data.IsFill(ix) ? fill : op(ix);
-   WriteVar(name, buf);
+   if(tdep)
+    WriteVar(name, tind, buf);
+   else
+    WriteVar(name, buf);
   }
   else
   {
@@ -427,7 +435,10 @@ class NCFileW: public NCFileWBase
    float        buf[c[0] * c[1]];
    for(size_t iy = 0; iy < c[0]; iy++)
     for(size_t ix = 0; ix < c[1]; ix++) buf[iy * c[1] + ix] = data.IsFill(ix, iy) ? fill : op(ix, iy);
-   WriteVar(name, buf);
+   if(tdep)
+    WriteVar(name, tind, buf);
+   else
+    WriteVar(name, buf);
   }
 
   if(!*this) return "Can't write variable " + name + ": " + ErrMessage();
@@ -435,14 +446,24 @@ class NCFileW: public NCFileWBase
   return "";
  }
 
- template<class D> MString WriteVariable(const D& data, const MString& name)
+ template<class D> MString WriteVariable(const D& data, const MString& name, size_t tind)
  {
   if constexpr(Is1DType(DetType<D>()))
-   return WriteVariable(data, name, [&data = std::as_const(data)](size_t i) { return data(i); });
+   return WriteVariable(
+       data, name, [&data = std::as_const(data)](size_t i) { return data(i); }, tind);
   else
-   return WriteVariable(data, name, [&data = std::as_const(data)](size_t i, size_t j) { return data(i, j); });
+   return WriteVariable(
+       data, name, [&data = std::as_const(data)](size_t i, size_t j) { return data(i, j); }, tind);
  }
 
+ template<class D, class Op>
+ requires(std::is_invocable_v<Op, size_t> || std::is_invocable_v<Op, size_t, size_t>)
+ MString WriteVariable(const D& data, const MString& name, Op op)
+ {
+  return WriteVariable(data, name, op, 0);
+ }
+ template<class D> MString WriteVariable(const D& data, const MString& name) { return WriteVariable(data, name, 0); }
+
  template<class D> MString WriteGrid(const D& data)
  {
   static constexpr auto dtype = DetType<D>();

src/ncfilew.cpp

@@ -87,16 +87,48 @@ MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, int compre
  return "";
 }
 
+MString NCFileW::AddTimeData(const TimeData& tdata, bool tisindex)
+{
+ tdep = tisindex;
+
+ AddDim("time", tdata.steps.size());
+ AddVar("time", Type2NCType<decltype(tdata.steps)::value_type>, "time");
+ SetComp("time", compress);
+ AddAtt("time", "standard_name", "time");
+ AddAtt("time", "long_name", "time");
+ AddAtt("time", "units", tdata.UnitName());
+ WriteVar("time", tdata.steps.data());
+
+ if(!*this) return MString("Can't add time data to the netcdf file: ") + ErrMessage();
+
+ return "";
+}
+
 MString NCFileW::AddVariable(const MString& name, const MString& stname, const MString& lname, const MString& units, const MString& comment)
 {
  if(type == UNKNOWN) return "File not opened";
 
  if(Is1DType(type))
-  AddVar(name, NC_FLOAT, "i");
+ {
+  if(tdep)
+   AddVar(name, NC_FLOAT, "time", "i");
+  else
+   AddVar(name, NC_FLOAT, "i");
+ }
  else if(IsGeoType(type))
-  AddVar(name, NC_FLOAT, "latitude", "longitude");
+ {
+  if(tdep)
+   AddVar(name, NC_FLOAT, "time", "latitude", "longitude");
+  else
+   AddVar(name, NC_FLOAT, "latitude", "longitude");
+ }
  else
-  AddVar(name, NC_FLOAT, "y", "x");
+ {
+  if(tdep)
+   AddVar(name, NC_FLOAT, "time", "y", "x");
+  else
+   AddVar(name, NC_FLOAT, "y", "x");
+ }
 
  SetComp(name, compress);
  if(stname.Exist()) AddAtt(name, "standard_name", stname);