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

Also, one-index version of Read() function.

actions/actiontsc.h

@@ -20,6 +20,7 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
  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;
@@ -40,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))};
  {
@@ -72,59 +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, args.at("_cmdline"), compress);
-  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

@@ -1,13 +1,14 @@
 #define MICHLIB_NOSOURCE
 #include "actionuv.h"
 
-void UVMethods::StPoints::WriteBinBile(const MString& name) const
+void UVMethods::StPoints::WriteBinBile(const MString& name, const michlib_internal::ParameterListEx& pars) const
 {
  BFileW stp;
  stp.Create(name, 3);
  stp.SetColumnName(1, lonlat ? "Longitude" : "x");
  stp.SetColumnName(2, lonlat ? "Latitude" : "x");
  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)");
+ stp.SetParameters(pars);
 
  for(size_t i = 0; i < N(); i++)
  {
@@ -19,9 +20,9 @@ void UVMethods::StPoints::WriteBinBile(const MString& name) const
  stp.Close();
 }
 
-MString UVMethods::StPoints::WriteNcFile(const MString& name, 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";
@@ -30,9 +31,26 @@ MString UVMethods::StPoints::WriteNcFile(const MString& name, const MString& his
  if(!nc) return "Can't create netcdf file " + name + ": " + nc.ErrMessage();
 
  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);
 
@@ -49,20 +67,38 @@ MString UVMethods::StPoints::WriteNcFile(const MString& name, const MString& his
   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 WriteBinBile(const MString& name) const;
+  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 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(); }
  };
 };
 
@@ -105,12 +117,18 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
  michlib_internal::ParameterListEx pars;
  pars.UsePrefix("");
  pars.SetParameter("source", args.at("source"));
+ pars.SetParameter("history", args.at("_cmdline"));
 
  bool mode = args.contains("geostrophic");
 
  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>)
  {
@@ -129,157 +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, args.at("_cmdline"), compress);
-   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);
-
-   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, args.at("_cmdline"), compress);
-   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);
-  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, 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/actiondep.h

@@ -5,6 +5,8 @@
 #include "mregex.h"
 #include "traits.h"
 #include "uvdata.h"
+#include <array>
+#include <utility>
 
 using michlib::MDateTime;
 
@@ -21,6 +23,41 @@ using michlib::MDateTime;
  };
 #endif
 
+struct TimeData
+{
+ static constexpr auto                                            stepunits = std::to_array<uint>({3600 * 24, 3600, 60, 1});
+ static constexpr std::array<const std::string, stepunits.size()> stepnames{"days", "hours", "minutes", "seconds"};
+
+ MDateTime                  refdate;
+ size_t                     unitind;
+ std::vector<michlib::int4> steps;
+
+ template<class D> TimeData(const D& data, const TIndex& tindexes): refdate(), unitind(stepunits.size()), steps(tindexes.size())
+ {
+  if(steps.size() == 0) return;
+  refdate = data.Time(tindexes[0]);
+  unitind = 0;
+
+  for(size_t i = 0; i < steps.size(); i++)
+  {
+   auto delta = data.Time(tindexes[i]) - refdate;
+   while(delta % stepunits[unitind] != 0) unitind++;
+  }
+
+  for(size_t i = 0; i < steps.size(); i++)
+  {
+   auto delta = data.Time(tindexes[i]) - refdate;
+   steps[i]   = michlib::int_cast<decltype(steps)::value_type>(delta / stepunits[unitind]);
+  }
+ }
+
+ MString UnitName() const
+ {
+  MString out = stepnames[unitind].c_str();
+  return out + " since " + refdate.ToString();
+ }
+};
+
 template<class D> size_t GetTIndex(const D& data, const MDateTime& t)
 {
  size_t nt = data.NTimes();
@@ -45,21 +82,27 @@ template<class D> std::pair<TIndex, MString> GetTIndexes(const D& data, const CL
   MString   regex = args.at("time");
   MDateTime time;
 
-  if(time.FromString(regex)) return {TIndex(1, GetTIndex(data, time)), ""};             // Time, not regex
-  if(regex == "BEGIN" || regex == "BEG" || regex == "FIRST") return {TIndex(1, 0), ""}; // First time
-  if(regex == "END" || regex == "LAST") return {TIndex(1, data.NTimes() - 1), ""};      // Last time
-
-  michlib::RegExpSimple reg(regex.Buf());
-  if(reg.Compile() != 0) return {tindexes, "Bad regular expression: " + regex};
-
-  for(size_t i = 0; i < data.NTimes(); i++)
+  if(time.FromString(regex)) // Time, not regex
+   tindexes.push_back(GetTIndex(data, time));
+  else if(regex == "BEGIN" || regex == "BEG" || regex == "FIRST") // First time
+   tindexes.push_back(0);
+  else if(regex == "END" || regex == "LAST") // Last time
+   tindexes.push_back(data.NTimes() - 1);
+  else // Regular expression
   {
-   MString date = data.Time(i).ToString();
-   if(reg.Match(date.Buf())) tindexes.push_back(i);
+   michlib::RegExpSimple reg(regex.Buf());
+   if(reg.Compile() != 0) return {tindexes, "Bad regular expression: " + regex};
+
+   for(size_t i = 0; i < data.NTimes(); i++)
+   {
+    MString date = data.Time(i).ToString();
+    if(reg.Match(date.Buf())) tindexes.push_back(i);
+   }
   }
+
   if(tindexes.size() == 0) return {tindexes, "There are no times matching the regular expression: " + regex};
   if(tindexes.size() == 1)
-   pars.SetParameter("time", data.Time(tindexes[0]).ToString());
+   pars.SetParameter("time", data.Time(tindexes[0]).ToTString());
   else
    pars.SetParameter("timeregex", args.at("time"));
  }
@@ -79,53 +122,46 @@ template<class D> std::pair<TIndex, MString> GetTIndexes(const D& data, const CL
   if(beg > data.Time(nt - 1)) return {tindexes, "Begin time " + b.ToTString() + " is greater then end time in the dataset " + data.Time(nt - 1).ToTString()};
   if(end < data.Time(0)) return {tindexes, "End time " + e.ToTString() + " is lesser then begin time in the dataset " + data.Time(0).ToTString()};
 
-  size_t ib = 0, ie = nt - 1;
   for(size_t i = 0; i < nt; i++)
-   if(data.Time(i) >= beg)
-   {
-    ib = i;
-    break;
-   }
-
-  for(size_t i = nt; i != 0; i--)
-   if(data.Time(i - 1) <= end)
-   {
-    ie = i - 1;
-    break;
-   }
-
-  tindexes.resize(ie - ib + 1);
-  for(size_t i = 0; i < ie - ib + 1; i++) tindexes[i] = i + ib;
+   if(data.Time(i) >= beg && data.Time(i) <= end) tindexes.push_back(i);
 
   if(tindexes.size() == 0) return {tindexes, "There are no times between " + b.ToString() + " and " + e.ToString()};
-  pars.SetParameter("timeb", b.ToString());
-  pars.SetParameter("timee", e.ToString());
+  pars.SetParameter("timeb", b.ToTString());
+  pars.SetParameter("timee", e.ToTString());
  }
 
+ std::ranges::sort(tindexes, [&data = std::as_const(data)](size_t a, size_t b) { return data.Time(a) < data.Time(b); });
+
  return {tindexes, ""};
 }
 
+template<class D> std::vector<ReadType<D>> Read(const D& data, const std::vector<MString>& vnames, const BaseParameters* p, size_t ind)
+{
+ using RT = ReadType<D>;
+ std::vector<RT> out;
+ michlib::message("Time: " + data.Time(ind).ToTString());
+ std::map<MString, RT> cache;
+ for(const auto& vname: vnames)
+ {
+  bool res;
+  if constexpr(ReadPSupported<D>)
+   res = data.Read(vname, cache, p, ind);
+  else if constexpr(ReadSupported<D>)
+   res = data.Read(vname, cache, ind);
+  if(!res) return out;
+ }
+ for(size_t i = 0; i < vnames.size(); i++) out.emplace_back(std::move(cache[vnames[i]]));
+
+ return out;
+}
+
 template<class D> std::vector<ReadType<D>> Read(const D& data, const std::vector<MString>& vnames, const BaseParameters* p, const TIndex& tindex)
 {
  using RT = ReadType<D>;
  size_t          ind;
  std::vector<RT> out;
  if(tindex.size() == 1)
- {
-  ind = tindex[0];
-  michlib::message("Time: " + data.Time(ind).ToTString());
-  std::map<MString, RT> cache;
-  for(const auto& vname: vnames)
-  {
-   bool res;
-   if constexpr(ReadPSupported<D>)
-    res = data.Read(vname, cache, p, ind);
-   else if constexpr(ReadSupported<D>)
-    res = data.Read(vname, cache, ind);
-   if(!res) return out;
-  }
-  for(size_t i = 0; i < vnames.size(); i++) out.emplace_back(std::move(cache[vnames[i]]));
- }
+  return Read(data, vnames, p, tindex[0]);
  else
  {
   std::vector<Averager<RT>> aver(vnames.size());

include/layereddata.h

@@ -49,14 +49,15 @@ class LayeredData: public NCFuncs
   }
 
   const NCFileA* operator->() const { return &nc; }
-  explicit operator bool() const { return nc; }
-  MDateTime Begin() const { return times.front(); }
-  MDateTime End() const { return times.back(); }
+  explicit       operator bool() const { return nc; }
+  MDateTime      Begin() const { return times.front(); }
+  MDateTime      End() const { return times.back(); }
 
   const NCFileA& Get() const { return nc; }
 
   const std::vector<MDateTime>& Times() const { return times; }
-  size_t                        Index(MDateTime tm) const
+
+  size_t Index(MDateTime tm) const
   {
    if(tm < Begin() || tm > End()) return 0;
    size_t b = 0, e = times.size() - 1;
@@ -73,6 +74,8 @@ class LayeredData: public NCFuncs
    }
    return 0;
   }
+
+  const MString& Url() const { return url; }
  };
  std::vector<NC>        nc;
  std::vector<real>      depths;

include/ncfilew.h

@@ -1,5 +1,6 @@
 #pragma once
 #include "MString.h"
+#include "actiondep.h"
 #include "traits.h"
 #include <algorithm>
 #include <netcdf.h>
@@ -23,26 +24,31 @@ class NCFileWBase
  {
   static constexpr nc_type nc = NC_UBYTE;
   static int               put_var(int nc, int vid, const michlib::uint1* data) { return nc_put_var_ubyte(nc, vid, data); }
+  static int put_vara(int nc, int vid, const size_t* startp, const size_t* countp, const michlib::uint1* data) { return nc_put_vara_ubyte(nc, vid, startp, countp, data); }
  };
  template<class Dummy> struct NCTypeD<michlib::int2, Dummy>
  {
   static constexpr nc_type nc = NC_SHORT;
   static int               put_var(int nc, int vid, const michlib::int2* data) { return nc_put_var_short(nc, vid, data); }
+  static int put_vara(int nc, int vid, const size_t* startp, const size_t* countp, const michlib::int2* data) { return nc_put_vara_short(nc, vid, startp, countp, data); }
  };
  template<class Dummy> struct NCTypeD<michlib::uint2, Dummy>
  {
   static constexpr nc_type nc = NC_USHORT;
   static int               put_var(int nc, int vid, const michlib::uint2* data) { return nc_put_var_ushort(nc, vid, data); }
+  static int put_vara(int nc, int vid, const size_t* startp, const size_t* countp, const michlib::uint2* data) { return nc_put_vara_ushort(nc, vid, startp, countp, data); }
  };
  template<class Dummy> struct NCTypeD<michlib::int4, Dummy>
  {
   static constexpr nc_type nc = NC_INT;
   static int               put_var(int nc, int vid, const michlib::int4* data) { return nc_put_var_int(nc, vid, data); }
+  static int put_vara(int nc, int vid, const size_t* startp, const size_t* countp, const michlib::int4* data) { return nc_put_vara_int(nc, vid, startp, countp, data); }
  };
  template<class Dummy> struct NCTypeD<michlib::uint4, Dummy>
  {
   static constexpr nc_type nc = NC_UINT;
   static int               put_var(int nc, int vid, const michlib::uint4* data) { return nc_put_var_uint(nc, vid, data); }
+  static int put_vara(int nc, int vid, const size_t* startp, const size_t* countp, const michlib::uint4* data) { return nc_put_vara_uint(nc, vid, startp, countp, data); }
  };
  template<class Dummy> struct NCTypeD<michlib::int8, Dummy>
  {
@@ -56,11 +62,13 @@ class NCFileWBase
  {
   static constexpr nc_type nc = NC_FLOAT;
   static int               put_var(int nc, int vid, const float* data) { return nc_put_var_float(nc, vid, data); }
+  static int               put_vara(int nc, int vid, const size_t* startp, const size_t* countp, const float* data) { return nc_put_vara_float(nc, vid, startp, countp, data); }
  };
  template<class Dummy> struct NCTypeD<double, Dummy>
  {
   static constexpr nc_type nc = NC_DOUBLE;
   static int               put_var(int nc, int vid, const double* data) { return nc_put_var_double(nc, vid, data); }
+  static int               put_vara(int nc, int vid, const size_t* startp, const size_t* countp, const double* data) { return nc_put_vara_double(nc, vid, startp, countp, data); }
  };
  template<class Dummy> struct NCTypeD<MString, Dummy>
  {
@@ -84,7 +92,7 @@ class NCFileWBase
   public:
   constexpr Error(): err(NC_NOERR) {}
   constexpr Error(int n_err): err(n_err) {}
-  explicit operator bool() const { return err == NC_NOERR; }
+  explicit    operator bool() const { return err == NC_NOERR; }
   const char* ErrMessage() const { return nc_strerror(err); }
   int         Err() const { return err; }
   bool        IsErr() const { return err != NC_NOERR; }
@@ -100,7 +108,7 @@ class NCFileWBase
  const char* ErrMessage() const { return err.ErrMessage(); }
 
  explicit operator bool() const { return Ok(); }
- bool Ok() const { return opened && !err.IsErr(); }
+ bool     Ok() const { return opened && !err.IsErr(); }
  operator int() const { return ncid; }
 
  void Close()
@@ -124,11 +132,91 @@ class NCFileWBase
   err.Reset(nc_def_dim(ncid, name.Buf(), len, &dimid));
  }
 
+ MString AddAtts(const michlib_internal::ParameterListEx& pars)
+ {
+  for(const auto& p: pars.GetParameterList())
+  {
+   MString ret;
+   if(p.first.Prefix() != "") continue; // Only parameters with empty prefix considered for now
+   switch(p.second.valtype)
+   {
+   case michlib_internal::Parameter::UINT1:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.u1);
+    break;
+   }
+   case michlib_internal::Parameter::UINT2:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.u2);
+    break;
+   }
+   case michlib_internal::Parameter::UINT4:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.u4);
+    break;
+   }
+   case michlib_internal::Parameter::UINT8:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.u8);
+    break;
+   }
+   case michlib_internal::Parameter::INT1:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.i1);
+    break;
+   }
+   case michlib_internal::Parameter::INT2:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.i2);
+    break;
+   }
+   case michlib_internal::Parameter::INT4:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.i4);
+    break;
+   }
+   case michlib_internal::Parameter::INT8:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.i8);
+    break;
+   }
+   case michlib_internal::Parameter::FLOAT:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.r4);
+    break;
+   }
+   case michlib_internal::Parameter::DOUBLE:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.r8);
+    break;
+   }
+   case michlib_internal::Parameter::LDOUBLE:
+   {
+    ret = AddAtt(p.first.Name(), static_cast<double>(p.second.value.r10));
+    break;
+   }
+   case michlib_internal::Parameter::BOOL:
+   {
+    ret = AddAtt(p.first.Name(), p.second.value.b ? 1 : 0);
+    break;
+   }
+   case michlib_internal::Parameter::STRING:
+   {
+    ret = AddAtt(p.first.Name(), p.second.svalue);
+    break;
+   }
+   case michlib_internal::Parameter::INVALID: break; // Silently skip
+   }
+   if(ret.Exist()) return ret;
+  }
+  return "";
+ }
+
  template<class T>
  requires(Type2NCType<T> != NC_NAT)
- void AddAtt(const MString& name, const T& val)
+ MString AddAtt(const MString& name, const T& val)
  {
-  if(err.IsErr()) return;
+  if(err.IsErr()) return "Can't write attribute " + name + " due to previous errors";
   static constexpr auto nct    = Type2NCType<T>;
   static constexpr bool ischar = std::is_same_v<T, char*>;
 
@@ -138,19 +226,21 @@ class NCFileWBase
    err.Reset(nc_put_att_text(ncid, NC_GLOBAL, name.Buf(), strlen(val) + 1, val));
   else
    err.Reset(nc_put_att_text(ncid, NC_GLOBAL, name.Buf(), val.Len() + 1, val.Buf()));
+
+  return err.IsErr() ? "Can't write attribute " + name + ": " + ErrMessage() : "";
  }
 
  template<class T>
  requires(Type2NCType<T> != NC_NAT)
- void AddAtt(const MString& vname, const MString& name, const T& val)
+ MString AddAtt(const MString& vname, const MString& name, const T& val)
  {
-  if(err.IsErr()) return;
+  if(err.IsErr()) return "Can't write attribute " + name + " to the variable " + vname + " due to previous errors";
   static constexpr auto nct    = Type2NCType<T>;
   static constexpr bool ischar = std::is_same_v<T, char*>;
   int                   varid;
 
   err.Reset(nc_inq_varid(ncid, vname.Buf(), &varid));
-  if(err.IsErr()) return;
+  if(err.IsErr()) return "Can't find variable " + vname + ": " + ErrMessage();
 
   if constexpr(nct != NC_CHAR)
    err.Reset(nc_put_att(ncid, varid, name.Buf(), nct, 1, &val));
@@ -158,16 +248,19 @@ class NCFileWBase
    err.Reset(nc_put_att_text(ncid, varid, name.Buf(), strlen(val) + 1, val));
   else
    err.Reset(nc_put_att_text(ncid, varid, name.Buf(), val.Len() + 1, val.Buf()));
+
+  return err.IsErr() ? "Can't write attribute " + name + " to the variable " + vname + ": " + ErrMessage() : "";
  }
 
- void AddAtt(const MString& vname, const MString& name, const char* val)
+ MString AddAtt(const MString& vname, const MString& name, const char* val)
  {
-  if(err.IsErr()) return;
+  if(err.IsErr()) return "Can't write attribute " + name + " to the variable " + vname + " due to previous errors";
   int varid;
 
   err.Reset(nc_inq_varid(ncid, vname.Buf(), &varid));
-  if(err.IsErr()) return;
+  if(err.IsErr()) return "Can't find variable " + vname + ": " + ErrMessage();
   err.Reset(nc_put_att_text(ncid, varid, name.Buf(), strlen(val) + 1, val));
+  return err.IsErr() ? "Can't write attribute " + name + " to the variable " + vname + ": " + ErrMessage() : "";
  }
 
  template<class... T> void AddVar(const MString& vname, nc_type type, const T&... dnames) { AddVar(vname, type, std::vector<MString>({dnames...})); }
@@ -202,6 +295,40 @@ class NCFileWBase
   err.Reset(NCTypeD<T, void>::put_var(ncid, varid, data));
  }
 
+ template<class T>
+ requires requires(int nc, int vid, const size_t* start, const size_t* count, const T* d) {
+  {
+   NCTypeD<T, void>::put_vara(nc, vid, start, count, d)
+  } -> std::same_as<int>;
+ }
+ void WriteVar(const MString& vname, size_t ind, const T* data)
+ {
+  if(err.IsErr()) return;
+
+  int varid;
+  err.Reset(nc_inq_varid(ncid, vname.Buf(), &varid));
+  if(err.IsErr()) return;
+
+  int ndim;
+  err.Reset(nc_inq_var(ncid, varid, nullptr, nullptr, &ndim, nullptr, nullptr));
+  if(err.IsErr()) return;
+  std::vector<int> dimids(ndim);
+  err.Reset(nc_inq_var(ncid, varid, nullptr, nullptr, nullptr, dimids.data(), nullptr));
+  if(err.IsErr()) return;
+
+  std::vector<size_t> start(ndim), count(ndim);
+  start[0] = ind;
+  count[0] = 1;
+  for(size_t i = 1; i < dimids.size(); i++)
+  {
+   start[i] = 0;
+   err.Reset(nc_inq_dim(ncid, dimids[i], nullptr, &count[i]));
+   if(err.IsErr()) return;
+  }
+
+  err.Reset(NCTypeD<T, void>::put_vara(ncid, varid, start.data(), count.data(), data));
+ }
+
  private:
  // Members
  int   ncid;
@@ -240,6 +367,7 @@ class NCFileW: public NCFileWBase
 
  Type type = UNKNOWN;
  int  compress;
+ bool tdep = false;
 
  template<class D> static constexpr Type DetType()
  {
@@ -257,18 +385,18 @@ class NCFileW: public NCFileWBase
  static constexpr bool IsGeoType(Type t) { return t == GPSET || t == GGRID || t == GRGRID; }
  static constexpr bool Is1DType(Type t) { return t == PSET || t == GPSET; }
 
- MString CreateFile(Type stype, const MString& name, const MString& history, int compression, size_t nx, size_t ny);
+ MString CreateFile(Type stype, const MString& name, int compression, size_t nx, size_t ny);
 
  public:
  static constexpr auto Fill() { return fill; }
 
- template<class D> MString Create(const D& data, const MString& name, const MString& history, int compression)
+ template<class D> MString Create(const D& data, const MString& name, int compression)
  {
   if(type != UNKNOWN) return "File already created";
   if constexpr(Is1DType<D>())
-   return CreateFile(DetType<D>(), name, history, compression, data.N(), 0);
+   return CreateFile(DetType<D>(), name, compression, data.N(), 0);
   else
-   return CreateFile(DetType<D>(), name, history, compression, data.Nx(), data.Ny());
+   return CreateFile(DetType<D>(), name, compression, data.Nx(), data.Ny());
  }
 
  void Close()
@@ -277,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";
@@ -293,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
   {
@@ -301,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();
@@ -309,13 +446,23 @@ 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)
  {

src/layereddata.cpp

@@ -51,8 +51,9 @@ MString LayeredData::Open(const MString& dataset)
    nc.emplace_back(std::move(url));
    if(!nc.back())
    {
+    auto failedurl = nc.back().Url();
     nc.clear();
-    return "Can't connect to url " + url;
+    return "Can't connect to url " + failedurl;
    }
   }
   else

src/ncfilew.cpp

@@ -1,7 +1,7 @@
 #define MICHLIB_NOSOURCE
 #include "ncfilew.h"
 
-MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, const MString& history, int compression, size_t nx, size_t ny)
+MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, int compression, size_t nx, size_t ny)
 {
  if(stype == UNKNOWN) return "Can't determine file type";
 
@@ -12,7 +12,6 @@ MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, const MStr
  Open(name);
  if(!*this) return "Can't create netcdf file " + name + ": " + ErrMessage();
 
- AddAtt("history", history);
  AddAtt("node_offset", node_offset);
 
  switch(stype)
@@ -88,16 +87,48 @@ MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, const MStr
  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);