Added support for netcdf output for the uv action

actions/actionuv.cpp

@@ -0,0 +1,88 @@
+#define MICHLIB_NOSOURCE
+#include "actionuv.h"
+
+void UVMethods::StPoints::WriteBinBile(const MString& name) 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)");
+
+ for(size_t i = 0; i < N(); i++)
+ {
+  stp.Write(x[i]);
+  stp.Write(y[i]);
+  stp.Write(michlib::int_cast<size_t>(t[i]));
+ }
+ stp.Finalize();
+ stp.Close();
+}
+
+MString UVMethods::StPoints::WriteNcFile(const MString& name, const MString& history, int comp) const
+{
+ int     ret;
+ int     ncid;
+ int     dimid;
+ int     xid, yid, tid;
+ MString text;
+
+ ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &ncid);
+ if(ret != NC_NOERR) return "Can't create netcdf file: " + name;
+
+ ret = nc_put_att_text(ncid, NC_GLOBAL, "history", history.Len() + 1, history.Buf());
+ if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
+
+ ret = nc_def_dim(ncid, "i", N(), &dimid);
+ if(ret != NC_NOERR) return "Can't create dimension in the netcdf file";
+
+ ret = nc_def_var(ncid, lonlat ? "longitude" : "x", NC_FLOAT, 1, &dimid, &xid);
+ if(ret != NC_NOERR) return "Can't create x variable in the netcdf file";
+ ret = nc_def_var(ncid, lonlat ? "latitude" : "y", NC_FLOAT, 1, &dimid, &yid);
+ if(ret != NC_NOERR) return "Can't create y variable in the netcdf file";
+
+ if(lonlat)
+ {
+  text = "longitude";
+  ret  = nc_put_att_text(ncid, xid, "standard_name", text.Len() + 1, text.Buf());
+  if(ret != NC_NOERR) return "Can't write standard_name attribute of longitude variable in the netcdf file";
+  text = "latitude";
+  ret  = nc_put_att_text(ncid, yid, "standard_name", text.Len() + 1, text.Buf());
+  if(ret != NC_NOERR) return "Can't write standard_name attribute of latitude variable in the netcdf file";
+ }
+
+ ret = nc_def_var_deflate(ncid, xid, 1, 1, comp);
+ if(ret != NC_NOERR) return "Can't set deflate parameters for x variable in the netcdf file";
+ ret = nc_def_var_deflate(ncid, yid, 1, 1, comp);
+ if(ret != NC_NOERR) return "Can't set deflate parameters for y variable in the netcdf file";
+
+ ret = nc_def_var(ncid, "type", NC_UBYTE, 1, &dimid, &tid);
+ if(ret != NC_NOERR) return "Can't create type variable in the netcdf file";
+ ret = nc_def_var_deflate(ncid, tid, 1, 1, comp);
+ if(ret != NC_NOERR) return "Can't set deflate parameters for type variable in the netcdf file";
+ text = "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";
+ ret  = nc_put_att_text(ncid, tid, "long_name", text.Len() + 1, text.Buf());
+ if(ret != NC_NOERR) return "Can't write long_name attribute of type variable in the netcdf file";
+
+ ret = nc_enddef(ncid);
+ if(ret != NC_NOERR) return "Can't finish definition of the netcdf file";
+
+ const size_t i = 0, c = N();
+ float        buf[c];
+
+ for(size_t ix = 0; ix < c; ix++) buf[ix] = x[ix];
+ ret = nc_put_vara_float(ncid, xid, &i, &c, buf);
+ if(ret != NC_NOERR) return "Can't write x variable in the netcdf file";
+
+ for(size_t ix = 0; ix < c; ix++) buf[ix] = y[ix];
+ ret = nc_put_vara_float(ncid, yid, &i, &c, buf);
+ if(ret != NC_NOERR) return "Can't write y variable in the netcdf file";
+
+ ret = nc_put_vara_ubyte(ncid, tid, &i, &c, t.data());
+ if(ret != NC_NOERR) return "Can't write type variable in the netcdf file";
+
+ ret = nc_close(ncid);
+ if(ret != NC_NOERR) return "Can't close netcdf file";
+
+ return "";
+}

actions/actionuv.h

@@ -1,11 +1,98 @@
 #pragma once
 #include "BFileW.h"
 #include "actiondep.h"
-#include <memory>
+#include "ncfilew.h"
+#include <utility>
 
 using michlib::BFileW;
 
-ADD_ACTION(UV, uv, ReadPSupported<Source> || ReadSupported<Source>);
+class UVMethods
+{
+ protected:
+ template<class D> class SparserHolder
+ {
+  protected:
+  const D& data;
+  size_t   shiftx, shifty, skipx, skipy;
+  SparserHolder(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): data(d), shiftx(shx), shifty(shy), skipx(sx), skipy(sy) {}
+
+  size_t I2Ix(size_t i) const { return i * skipx + shiftx; }
+  size_t I2Iy(size_t i) const { return i * skipy + shifty; }
+
+  public:
+  size_t Nx() const
+  {
+   size_t n = data.Nx() - shiftx;
+   return n / skipx + ((n % skipx > 0) ? 1 : 0);
+  }
+  size_t Ny() const
+  {
+   size_t n = data.Ny() - shifty;
+   return n / skipy + ((n % skipy > 0) ? 1 : 0);
+  }
+
+  static constexpr auto Fillval() { return D::Fillval(); }
+
+  real Lon(size_t ix, size_t iy) const { return data.Lon(I2Ix(ix), I2Iy(iy)); }
+  real Lat(size_t ix, size_t iy) const { return data.Lat(I2Ix(ix), I2Iy(iy)); }
+  real U(size_t ix, size_t iy) const { return data.U(I2Ix(ix), I2Iy(iy)); }
+  real V(size_t ix, size_t iy) const { return data.V(I2Ix(ix), I2Iy(iy)); }
+  bool IsFill(size_t ix, size_t iy) const { return data.IsFill(I2Ix(ix), I2Iy(iy)); }
+ };
+
+ template<class D, bool b> class Sparser2DGeoRectArray: public SparserHolder<D>
+ {
+  protected:
+  Sparser2DGeoRectArray(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): SparserHolder<D>(d, shx, shy, sx, sy) {}
+ };
+ template<class D> class Sparser2DGeoRectArray<D, true>: public SparserHolder<D>
+ {
+  protected:
+  using SparserHolder<D>::data;
+  using SparserHolder<D>::I2Ix;
+  using SparserHolder<D>::I2Iy;
+
+  Sparser2DGeoRectArray(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): SparserHolder<D>(d, shx, shy, sx, sy) {}
+
+  public:
+  real Ix2Lon(size_t ix) const { return data.Ix2Lon(I2Ix(ix)); }
+  real Iy2Lat(size_t iy) const { return data.Iy2Lat(I2Iy(iy)); }
+ };
+
+ template<class D> class Sparser: public Sparser2DGeoRectArray<D, ReadIs2DGeoRectArray<D>>
+ {
+  public:
+  Sparser(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): Sparser2DGeoRectArray<D, ReadIs2DGeoRectArray<D>>(d, shx, shy, sx, sy) {}
+ };
+
+ class StPoints
+ {
+  std::vector<real>           x, y;
+  std::vector<michlib::uint1> t;
+  bool                        lonlat;
+
+  size_t N() const { return t.size(); }
+
+  public:
+  StPoints(bool ll): lonlat(ll) {}
+
+  void Add(const std::vector<UVDataStorage::StPoint>& points)
+  {
+   for(const auto& p: points)
+   {
+    x.push_back(p.x);
+    y.push_back(p.y);
+    t.push_back(michlib::int_cast<michlib::uint1>(p.type));
+   }
+  }
+
+  void WriteBinBile(const MString& name) const;
+
+  MString WriteNcFile(const MString& name, const MString& history, int comp) const;
+ };
+};
+
+ADD_ACTION(UV, uv, ReadPSupported<Source> || ReadSupported<Source>, UVMethods);
 
 template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
 {
@@ -45,20 +132,27 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
 
  // Main file
  MString name   = args.contains("out") ? args.at("out") : "";
+ MString outfmt = args.contains("outformat") ? args.at("outformat") : (GetExt(name) == "nc" ? "nc" : "bin");
+
+ MString u = data.Unit().Exist() ? data.Unit() : "unknown", d = data.DUnit().Exist() ? data.DUnit() : "unknown";
+
  if(name.Exist())
+ {
+  if(outfmt == "bin")
   {
    BFileW fw;
    fw.Create(name, 9);
    fw.SetColumnName(1, "Longitude");
    fw.SetColumnName(2, "Latitude");
-  fw.SetColumnName(3, "u, cm/s");
+   fw.SetColumnName(3, "u, " + u);
-  fw.SetColumnName(4, "v, cm/s");
+   fw.SetColumnName(4, "v, " + u);
-  fw.SetColumnName(5, "Div, (cm/s)/km");
+   fw.SetColumnName(5, "Div, (" + u + ")/" + d);
-  fw.SetColumnName(6, "Rot, (cm/s)/km");
+   fw.SetColumnName(6, "Rot, (" + u + ")/" + d);
-  fw.SetColumnName(7, "Okubo-Weiss parameter, (cm^2/s^2)/km^2");
+   fw.SetColumnName(7, "Okubo-Weiss parameter, (" + u + ")2/" + d + "2");
-  fw.SetColumnName(8, "Kinetic energy, cm^2/s^2");
+   fw.SetColumnName(8, "Kinetic energy, (" + u + ")2");
-  fw.SetColumnName(9, "Eddy kinetic energy, cm^2/s^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));
@@ -74,9 +168,44 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
    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;
+
+   fw.Close();
+  }
+  else
+   return "Unknown format: " + outfmt;
+ }
 
  // 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())
  {
   size_t shiftx = args.contains("shiftx") ? args.at("shiftx").ToInteger<size_t>() : 0;
@@ -84,48 +213,70 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
   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;
 
+  Sparser<decltype(data)> sdata(data, shiftx, shifty, skipx, skipy);
+  if(outfmt == "bin")
+  {
    BFileW vel;
    vel.Create(name, 4);
    vel.SetColumnName(1, "Longitude");
    vel.SetColumnName(2, "Latitude");
-  vel.SetColumnName(3, "u, cm/s");
+   vel.SetColumnName(3, "u, " + u);
-  vel.SetColumnName(4, "v, cm/s");
+   vel.SetColumnName(4, "v, " + u);
 
-  for(size_t ix = shiftx; ix < data.Nx(); ix += skipx)
+   for(size_t ix = 0; ix < sdata.Nx(); ix++)
-   for(size_t iy = shifty; iy < data.Ny(); iy += skipy)
+    for(size_t iy = 0; iy < sdata.Ny(); iy++)
     {
-    vel.Write(data.Lon(ix, iy));
+     vel.Write(sdata.Lon(ix, iy));
-    vel.Write(data.Lat(ix, iy));
+     vel.Write(sdata.Lat(ix, iy));
-    vel.Write(data.U(ix, iy) == data.Fillval() ? NAN : data.U(ix, iy));
+     vel.Write(sdata.U(ix, iy) == sdata.Fillval() ? NAN : sdata.U(ix, iy));
-    vel.Write(data.V(ix, iy) == data.Fillval() ? NAN : data.V(ix, iy));
+     vel.Write(sdata.V(ix, iy) == sdata.Fillval() ? NAN : sdata.V(ix, iy));
     }
    vel.Finalize();
    vel.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(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();
+  }
+  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())
  {
-  BFileW stp;
+  StPoints p(true);
-  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++)
+   for(size_t iy = 0; iy < data.Ny() - 1; iy++) p.Add(data.StablePoints(ix, iy));
-   {
+
-    auto points = data.StablePoints(ix, iy);
+  if(outfmt == "bin")
-    for(size_t i = 0; i < points.size(); i++)
+   p.WriteBinBile(name);
+  else if(outfmt == "nc" || outfmt == "netcdf")
   {
-     stp.Write(points[i].x);
+   int     compress = args.contains("compress") ? args.at("compress").ToInt() : 3;
-     stp.Write(points[i].y);
+   MString err      = p.WriteNcFile(name, args.at("_cmdline"), compress);
-     stp.Write(michlib::int_cast<size_t>(points[i].type));
+   if(err.Exist()) return err;
-    }
   }
-  stp.Finalize();
+  else
-  stp.Close();
+   return "Unknown format: " + outfmt;
  }
 
  return "";