Gradient calculation code (without filtering)

actions/actiongrad.cpp

@@ -0,0 +1,219 @@
+#define MICHLIB_NOSOURCE
+#include "actiongrad.h"
+
+MString GradMethods::NCFileW::Create(const MString& name, const MString& history, const std::vector<MString>& vnames, const std::vector<MString>& lnames,
+                                     const std::vector<GradMethods::DataType>& lons, const std::vector<GradMethods::DataType>& lats, int compress)
+{
+ int ret;
+
+ const float node_offset = 0.0;
+
+ auto nx = lons.size();
+ auto ny = lats.size();
+
+ int xid, yid;
+
+ MString text;
+
+ struct
+ {
+  union
+  {
+   struct
+   {
+    int ydimid, xdimid;
+   };
+   int dimid[2];
+  };
+ } dim;
+
+ struct Keeper
+ {
+  bool active;
+  int  ncid;
+  ~Keeper()
+  {
+   if(active) nc_close(ncid);
+  }
+  operator int() { return ncid; }
+ } newnc;
+
+ const auto fill = static_cast<GradMethods::DataType>(std::numeric_limits<GradMethods::Matrix::MDataType>::max());
+
+ // Creating
+ ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &newnc.ncid);
+ if(ret != NC_NOERR) return "Can't create netcdf file: " + name;
+ newnc.active = true;
+
+ // Global attributes
+ ret = nc_put_att_text(newnc, NC_GLOBAL, "history", history.Len() + 1, history.Buf());
+ if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
+
+ ret = nc_put_att(newnc, NC_GLOBAL, "node_offset", NC_FLOAT, 1, &node_offset);
+ if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
+
+ // Dimensions
+ ret = nc_def_dim(newnc, "longitude", nx, &dim.xdimid);
+ if(ret != NC_NOERR) return "Can't create x-dimension in the netcdf file";
+ ret = nc_def_dim(newnc, "latitude", ny, &dim.ydimid);
+ if(ret != NC_NOERR) return "Can't create y-dimension in the netcdf file";
+
+ // lon, lat variables
+ ret = nc_def_var(newnc, "longitude", NC_FLOAT, 1, &dim.xdimid, &xid);
+ if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file";
+ ret = nc_def_var(newnc, "latitude", NC_FLOAT, 1, &dim.ydimid, &yid);
+ if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file";
+
+ ret = nc_def_var_deflate(newnc, xid, 1, 1, compress);
+ if(ret != NC_NOERR) return "Can't set deflate parameters for longitude variable in the netcdf file";
+ ret = nc_def_var_deflate(newnc, yid, 1, 1, compress);
+ if(ret != NC_NOERR) return "Can't set deflate parameters for latitude variable in the netcdf file";
+
+ text = "longitude";
+ ret  = nc_put_att_text(newnc, 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(newnc, 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";
+
+ text = "Longitude";
+ ret  = nc_put_att_text(newnc, xid, "long_name", text.Len() + 1, text.Buf());
+ if(ret != NC_NOERR) return "Can't write long_name attribute of longitude variable in the netcdf file";
+ text = "Latitude";
+ ret  = nc_put_att_text(newnc, yid, "long_name", text.Len() + 1, text.Buf());
+ if(ret != NC_NOERR) return "Can't write long_name attribute of latitude variable in the netcdf file";
+
+ // Variables
+ for(size_t i = 0; i < vnames.size(); i++)
+ {
+  int vid;
+  ret = nc_def_var(newnc, vnames[i].Buf(), NC_FLOAT, 2, dim.dimid, &vid);
+  if(ret != NC_NOERR) return "Can't create " + vnames[i] + " variable in the netcdf file";
+
+  ret = nc_def_var_deflate(newnc, vid, 1, 1, compress);
+  if(ret != NC_NOERR) return "Can't set deflate parameters for " + vnames[i] + " variable in the netcdf file";
+
+  if(lnames[i].Exist()) ret = nc_put_att_text(newnc, vid, "long_name", lnames[i].Len() + 1, lnames[i].Buf());
+  if(ret != NC_NOERR) return "Can't write long_name attribute of " + vnames[i] + " variable in the netcdf file";
+
+  ret = nc_put_att_float(newnc, vid, "_FillValue", NC_FLOAT, 1, &fill);
+  if(ret != NC_NOERR) return "Can't write _FillValue attribute of " + vnames[i] + " variable in the netcdf file";
+ }
+
+ // End definitions
+ nc_enddef(newnc);
+
+ // Writing lon, lat
+ const size_t i = 0;
+
+ ret = nc_put_vara_float(newnc, xid, &i, &nx, lons.data());
+ if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file";
+ ret = nc_put_vara_float(newnc, yid, &i, &ny, lats.data());
+ if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file";
+
+ ncid         = newnc;
+ newnc.active = false;
+ return "";
+}
+
+MString GradMethods::NCFileW::WriteVar(const MString& name, const std::vector<GradMethods::Matrix::MDataType>& data)
+{
+ int    ret;
+ int    vid;
+ int    xid, yid;
+ size_t nx, ny;
+
+ ret = nc_inq_varid(ncid, name.Buf(), &vid);
+ if(ret != NC_NOERR) return "Can't find variable " + name + " in the netcdf file";
+
+ ret = nc_inq_dimid(ncid, "longitude", &xid);
+ if(ret != NC_NOERR) return "Can't find longitude dimension in the netcdf file";
+ ret = nc_inq_dimid(ncid, "latitude", &yid);
+ if(ret != NC_NOERR) return "Can't find latitude dimension in the netcdf file";
+
+ ret = nc_inq_dimlen(ncid, xid, &nx);
+ if(ret != NC_NOERR) return "Can't find longitude size in the netcdf file";
+ ret = nc_inq_dimlen(ncid, yid, &ny);
+ if(ret != NC_NOERR) return "Can't find latitude size in the netcdf file";
+
+ const size_t i[2] = {0, 0};
+ const size_t c[2] = {ny, nx};
+
+ GradMethods::DataType buf[nx * ny];
+ for(size_t i = 0; i < nx * ny; i++) buf[i] = static_cast<GradMethods::DataType>(data[i]);
+
+ ret = nc_put_vara_float(ncid, vid, i, c, buf);
+ if(ret != NC_NOERR) return "Can't write " + name + " variable in the netcdf file";
+ return "";
+}
+
+GradMethods::Matrix::Matrix(const std::vector<GradMethods::DataType>& in, size_t nx_, size_t ny_, struct GradMethods::MinMax minmax): nx(nx_), ny(ny_), data(nx_ * ny_)
+{
+ if(minmax.automin || minmax.automax)
+ {
+  DataType min = in[0];
+  DataType max = in[0];
+  for(size_t i = 1; i < in.size(); i++)
+   if(in[i] != minmax.fill)
+   {
+    min = std::min(min, in[i]);
+    max = std::max(max, in[i]);
+   }
+  if(minmax.automin) minmax.min = min;
+  if(minmax.automax) minmax.max = max;
+ }
+
+ if(minmax.log)
+ {
+  minmax.min = michlib_internal::RealType<sizeof(DataType)>::Log(minmax.min);
+  minmax.max = michlib_internal::RealType<sizeof(DataType)>::Log(minmax.max);
+ }
+
+ DataType a = (std::numeric_limits<MDataType>::max() - 1) / (minmax.max - minmax.min);
+ for(size_t i = 1; i < in.size(); i++)
+ {
+  DataType v = minmax.log ? michlib_internal::RealType<sizeof(DataType)>::Log(in[i]) : in[i];
+  if(in[i] == minmax.fill)
+   data[i] = std::numeric_limits<MDataType>::max();
+  else if(v <= minmax.min)
+   data[i] = 0;
+  else if(v >= minmax.max)
+   data[i] = std::numeric_limits<MDataType>::max() - 1;
+  else
+   data[i] = static_cast<MDataType>(michlib_internal::RealType<sizeof(DataType)>::Round(a * (v - minmax.min)));
+ }
+}
+
+void GradMethods::Matrix::Grad()
+{
+ std::vector<MDataType> out(data.size());
+ const auto             bad = std::numeric_limits<MDataType>::max();
+
+ for(size_t iy = 0; iy < ny; iy++)
+  for(size_t ix = 0; ix < nx; ix++)
+  {
+   if(iy < 1 || ix < 1 || iy > ny - 2 || ix > nx - 2)
+    out[iy * nx + ix] = bad;
+   else if(V(ix - 1, iy - 1) == bad || V(ix, iy - 1) == bad || V(ix + 1, iy - 1) == bad || V(ix - 1, iy) == bad || V(ix, iy) == bad || V(ix + 1, iy) == bad ||
+           V(ix - 1, iy + 1) == bad || V(ix, iy + 1) == bad || V(ix + 1, iy + 1) == bad)
+    out[iy * nx + ix] = bad;
+   else
+   {
+    using IT = michlib::int4;
+    // Possible but unlikely overflow
+    const IT m1 = -1;
+    const IT m2 = -2;
+    const IT p1 = 1;
+    const IT p2 = 2;
+
+    IT gx = m1 * V(ix - 1, iy + 1) + p1 * V(ix + 1, iy + 1) + m2 * V(ix - 1, iy) + p2 * V(ix + 1, iy) + m1 * V(ix - 1, iy - 1) + p1 * V(ix + 1, iy - 1);
+    IT gy = m1 * V(ix - 1, iy - 1) + p1 * V(ix - 1, iy + 1) + m2 * V(ix, iy - 1) + p2 * V(ix, iy + 1) + m1 * V(ix + 1, iy - 1) + p1 * V(ix + 1, iy + 1);
+
+    auto sq = static_cast<IT>(michlib::Round(michlib::Hypot(gx, gy)));
+    if(sq >= bad) sq = bad - 1;
+    out[iy * nx + ix] = static_cast<MDataType>(sq);
+   }
+  }
+
+ data = std::move(out);
+}

actions/actiongrad.h

@@ -0,0 +1,146 @@
+#pragma once
+#include "actiondep.h"
+#include "ncfilew.h"
+#include "ncfuncs.h"
+
+class GradMethods
+{
+ public:
+ using DataType = float;
+
+ static DataType ToNum(const MString& str) { return michlib_internal::RealType<sizeof(DataType)>::String2Real(str.Buf()); }
+
+ struct MinMax
+ {
+  bool     automin, automax, log;
+  DataType min, max;
+  DataType fill;
+ };
+
+ class Matrix;
+ class NCFileW;
+};
+
+class GradMethods::Matrix
+{
+ public:
+ using MDataType = michlib::uint2;
+
+ private:
+ size_t                 nx, ny;
+ std::vector<MDataType> data;
+
+ public:
+ Matrix(const std::vector<DataType>& in, size_t nx_, size_t ny_, struct MinMax minmax);
+
+ void Grad();
+
+ auto Nx() const { return nx; }
+ auto Ny() const { return ny; }
+
+ const auto& V(size_t ix, size_t iy) const { return data[iy * nx + ix]; }
+
+ auto& V(size_t ix, size_t iy) { return data[iy * nx + ix]; }
+
+ const auto& Data() const { return data; }
+};
+
+class GradMethods::NCFileW
+{
+ bool opened;
+ int  ncid;
+
+ public:
+ NCFileW(): opened(false) {}
+
+ MString Create(const MString& name, const MString& history, const std::vector<MString>& vnames, const std::vector<MString>& lnames, const std::vector<GradMethods::DataType>& lons,
+                const std::vector<GradMethods::DataType>& lats, int compress);
+ MString WriteVar(const MString& name, const std::vector<GradMethods::Matrix::MDataType>& data);
+
+ ~NCFileW() { Close(); }
+
+ void Close()
+ {
+  if(opened) nc_close(ncid);
+  opened = false;
+ }
+};
+
+template<class T>
+concept GradSupported = requires(T t, const MString& vname) {
+ {
+  t.ReadVar(vname)
+ } -> std::same_as<std::vector<GradMethods::DataType>>;
+};
+
+ADD_ACTION(GRAD, grad, GradSupported<Source>, GradMethods);
+
+template<class D> MString ActionGRAD::DoAction(const CLArgs& args, D& ds)
+{
+ auto resop = ds.Open(args);
+ if(resop.Exist()) return "Can't open source: " + resop;
+
+ MString name = args.contains("out") ? args.at("out") : "out.nc";
+
+ MString min = args.contains("min") ? args.at("min") : "auto";
+ MString max = args.contains("max") ? args.at("max") : "auto";
+
+ int compress = args.contains("compress") ? args.at("compress").ToInt() : 3;
+
+ std::vector<Matrix> data;
+
+ std::vector<MString> lnames;
+
+ // Read data
+ for(size_t i = 0; i < ds.NVar(); i++)
+ {
+  const MString& name  = ds.VarNames()[i];
+  const MString& lname = ds.LongNames()[i];
+  bool           hmin  = args.contains(name + "_min");
+  bool           hmax  = args.contains(name + "_max");
+  struct MinMax  minmax;
+
+  minmax.log = args.contains(name + "_log");
+
+  if(hmin)
+  {
+   MString vmin   = args.at(name + "_min");
+   minmax.automin = (vmin == "auto");
+   minmax.min     = ToNum(vmin);
+  }
+  else
+  {
+   minmax.automin = (min == "auto");
+   minmax.min     = ToNum(min);
+  }
+
+  if(hmax)
+  {
+   MString vmax   = args.at(name + "_max");
+   minmax.automax = (vmax == "auto");
+   minmax.max     = ToNum(vmax);
+  }
+  else
+  {
+   minmax.automax = (max == "auto");
+   minmax.max     = ToNum(max);
+  }
+
+  minmax.fill = ds.FillVal(name);
+
+  data.emplace_back(ds.ReadVar(name), ds.Nx(), ds.Ny(), minmax);
+  lnames.emplace_back(lname + ", gradient");
+ }
+
+ NCFileW fw;
+ fw.Create(name, (ds.History().Exist() ? (ds.History() + "; ") : "") + args.at("_cmdline"), ds.VarNames(), lnames, ds.ReadLons(), ds.ReadLats(), compress);
+
+ for(size_t i = 0; i < ds.NVar(); i++)
+ {
+  const MString& name = ds.VarNames()[i];
+  data[i].Grad();
+  fw.WriteVar(name, data[i].Data());
+ }
+
+ return "";
+};

sources/TSCDATA.cpp

@@ -0,0 +1,129 @@
+#define MICHLIB_NOSOURCE
+#include "TSCDATA.h"
+
+MString TSCDATAData::Open(const CLArgs& args)
+{
+ if(!args.contains("dataset")) return "path to data not specified";
+ MString dataset = args.at("dataset");
+
+ michlib::NCFileA     newnc;
+ std::vector<MString> newvnames, newlnames;
+ MString              newhistory;
+
+ newnc.Reset(dataset);
+ if(!newnc) return "Can't open file " + dataset;
+
+ auto head = newnc.Header();
+ if(head.Dimensions().size() != 2) return "Unsupported number of dimensions";
+ if((head.Dimensions()[0].Name() != "longitude" || head.Dimensions()[1].Name() != "latitude") &&
+    (head.Dimensions()[1].Name() != "longitude" || head.Dimensions()[0].Name() != "latitude"))
+  return "Unsupported dimensions names";
+
+ if(head.Dimensions()[0].Name() == "longitude")
+ {
+  nx = head.Dimensions()[0].Len();
+  ny = head.Dimensions()[1].Len();
+ }
+ else
+ {
+  ny = head.Dimensions()[0].Len();
+  nx = head.Dimensions()[1].Len();
+ }
+
+ {
+  bool lonfound = false, latfound = false;
+  for(const auto& v: head.Variables())
+   if(v.Dimensions().size() == 1 && v.Type().Id() == NC_FLOAT)
+   {
+    lonfound = lonfound || (v.Dimensions()[0].Name() == "longitude" && v.Name() == "longitude");
+    latfound = latfound || (v.Dimensions()[0].Name() == "latitude" && v.Name() == "latitude");
+   }
+  if(!lonfound) return "Longitude not found";
+  if(!latfound) return "Latitude not found";
+ }
+
+ for(const auto& v: head.Variables())
+ {
+  if(v.Dimensions().size() != 2) continue;
+  if(v.Type().Id() != NC_FLOAT) continue;
+  if((v.Dimensions()[0].Name() != "longitude" || v.Dimensions()[1].Name() != "latitude") && (v.Dimensions()[1].Name() != "longitude" || v.Dimensions()[0].Name() != "latitude"))
+   continue;
+  newvnames.push_back(v.Name());
+  auto lname = newnc.A<MString>(v.Name(), "long_name");
+  newlnames.push_back(lname ? lname.Get() : "");
+ }
+ if(newvnames.size() == 0) return "No variables found";
+
+ {
+  auto his = newnc.A<MString>("history");
+  if(his) newhistory = his;
+ }
+
+ history = std::move(newhistory);
+ vnames  = std::move(newvnames);
+ lnames  = std::move(newlnames);
+ nc      = std::move(newnc);
+ return "";
+}
+
+MString TSCDATAData::Info() const
+{
+ if(!nc) return "";
+
+ MString out;
+
+ out += MString("Dimensions: ") + nx + " X " + ny + "\n";
+ out += MString("Variables: ");
+ for(size_t i = 0; i < vnames.size(); i++) out += ((i == 0) ? "" : ", ") + vnames[i];
+ out += "\n";
+ auto his = nc.A<MString>("history");
+ if(his) out += "Creator: " + his.Get() + "\n";
+ return out;
+}
+
+std::vector<TSCDATAData::DataType> TSCDATAData::ReadLons() const
+{
+ std::vector<DataType> out;
+ auto                  lons = nc.V<DataType>("longitude");
+ if(lons)
+ {
+  out.resize(lons.DimLen(0));
+  for(size_t i = 0; i < out.size(); i++) out[i] = lons(i);
+ }
+ return out;
+}
+
+std::vector<TSCDATAData::DataType> TSCDATAData::ReadLats() const
+{
+ std::vector<DataType> out;
+ auto                  lats = nc.V<DataType>("latitude");
+ if(lats)
+ {
+  out.resize(lats.DimLen(0));
+  for(size_t i = 0; i < out.size(); i++) out[i] = lats(i);
+ }
+ return out;
+}
+
+std::vector<TSCDATAData::DataType> TSCDATAData::ReadVar(const MString& name) const
+{
+ std::vector<DataType> out;
+ bool                  havevar = false;
+ for(size_t i = 0; i < vnames.size(); i++) havevar = havevar || (vnames[i] == name);
+ if(!havevar) return out;
+
+ auto var = nc.V<DataType>(name, "latitude", "longitude");
+ if(var)
+ {
+  out.resize(nx * ny);
+  for(size_t iy = 0; iy < ny; iy++)
+   for(size_t ix = 0; ix < nx; ix++) out[iy * nx + ix] = var(iy, ix);
+ }
+ return out;
+}
+
+TSCDATAData::DataType TSCDATAData::FillVal(const MString& name) const
+{
+ auto fill = nc.A<DataType>(name, "_FillValue");
+ return fill ? fill.Get() : 0.0;
+}

sources/TSCDATA.h

@@ -0,0 +1,35 @@
+#pragma once
+#include "DataAdapters/ncfilealt.h"
+#include "ParseArgs.h"
+#include "simple2ddata.h"
+
+class TSCDATAData
+{
+ michlib::NCFileA     nc;
+ std::vector<MString> vnames, lnames;
+ size_t               nx, ny;
+ MString              history;
+
+ using DataType = float;
+
+ public:
+ static constexpr const char* name = "TSCDATA";
+
+ MString Info() const;
+ MString Open(const CLArgs& args);
+
+ std::vector<DataType> ReadLons() const;
+ std::vector<DataType> ReadLats() const;
+ std::vector<DataType> ReadVar(const MString& name) const;
+
+ const auto& VarNames() const { return vnames; }
+ const auto& LongNames() const { return lnames; }
+
+ size_t Nx() const { return nx; }
+ size_t Ny() const { return ny; }
+ size_t NVar() const { return vnames.size(); }
+
+ const auto& History() const { return history; }
+
+ DataType FillVal(const MString& name) const;
+};