Calculatin gradients in the action tsc

actions/actiontsc.h

@@ -41,7 +41,11 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
   }
  }
 
- bool average = args.contains("average");
+ bool average  = args.contains("average");
+ bool gradient = args.contains("gradient");
+
+ if(gradient)
+  if constexpr(!ReadIs2DGeoArray<D>) return "Gradient calculation not supported for this source";
 
  int compress = 3;
  if(args.contains("compress")) compress = args.at("compress").ToInt();
@@ -108,13 +112,24 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
    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"));
+   for(size_t i = 0; i < data.size(); i++)
+    fw.SetColumnName(++ic, (gradient ? "gradient of " : "") + vnames[i] + ", " + (data[i].Unit().Exist() ? data[i].Unit() : "unknown") + (gradient ? "/km" : ""));
    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));
+    if(gradient)
+    {
+     if constexpr(ReadIs2DGeoArray<D>)
+      for(size_t i = 0; i < data.size(); i++)
+      {
+       auto val = data[i].Grad(r);
+       fw.Write(val == data[i].Fillval() ? NAN : val);
+      }
+    }
+    else
+     for(size_t i = 0; i < data.size(); i++) fw.Write(data[i].IsFill(r) ? NAN : data[i](r));
    }
    fw.Finalize();
    fw.Close();
@@ -127,10 +142,21 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
    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.AddVariable((gradient ? "G" : "") + vnames[i], data[i].StandartName(), data[i].LongName(), data[i].Unit() + (gradient ? "/km" : ""), 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())
+    {
+     if(gradient)
+     {
+      if constexpr(ReadIs2DGeoArray<D>)
+       err = average ? ncfw.WriteVariable(data[i], "G" + vnames[i], [&data = std::as_const(data[i])](size_t i, size_t j) { return data.Grad(i, j); })
+                     : ncfw.WriteVariable(data[i], "G" + vnames[i], [&data = std::as_const(data[i])](size_t i, size_t j) { return data.Grad(i, j); }, it);
+     }
+     else
+      err = average ? ncfw.WriteVariable(data[i], vnames[i]) : ncfw.WriteVariable(data[i], vnames[i], it);
+    }
    if(err.Exist()) return err;
   }
   else

include/simple2ddata.h

@@ -1,5 +1,8 @@
 #pragma once
 #include "basedata.h"
+#include "geohelpers.h"
+
+using michlib::M_PI;
 
 class Simple2DData: public BaseData
 {
@@ -7,6 +10,8 @@ class Simple2DData: public BaseData
  size_t nx = 0, ny = 0;
  real   xstep = 0.0, ystep = 0.0;
 
+ static real D(real lon1, real lat1, real lon2, real lat2) { return michlib::GCD(lon1 * M_PI / 180.0, lat1 * M_PI / 180.0, lon2 * M_PI / 180.0, lat2 * M_PI / 180.0) * 6371.0; }
+
  public:
  using BaseData::IsFill;
  using BaseData::V;
@@ -59,12 +64,66 @@ class Simple2DData: public BaseData
 
   return out;
  }
+
+ real dVdX(size_t ix, size_t iy) const
+ {
+  if(IsFill(ix, iy)) return Fillval();
+  if(ix == 0 || IsFill(ix - 1, iy))
+  {
+   if(IsFill(ix + 1, iy))
+    return Fillval();
+   else
+    return (V(ix + 1, iy) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy));
+  }
+  if(ix == nx - 1 || IsFill(ix + 1, iy))
+  {
+   if(IsFill(ix - 1, iy))
+    return Fillval();
+   else
+    return (V(ix, iy) - V(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy));
+  }
+  return 0.5 * ((V(ix + 1, iy) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
+                (V(ix, iy) - V(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
+ }
+ real dVdY(size_t ix, size_t iy) const
+ {
+  if(IsFill(ix, iy)) return Fillval();
+
+  if(iy == 0 || IsFill(ix, iy - 1))
+  {
+   if(IsFill(ix, iy + 1))
+    return Fillval();
+   else
+    return (V(ix, iy + 1) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1));
+  }
+  if(iy == ny - 1 || IsFill(ix, iy + 1))
+  {
+   if(IsFill(ix, iy - 1))
+    return Fillval();
+   else
+    return (V(ix, iy) - V(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy));
+  }
+  return 0.5 * ((V(ix, iy + 1) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
+                (V(ix, iy) - V(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
+ }
+
+ real Grad(size_t ix, size_t iy) const
+ {
+  real grx = dVdX(ix, iy);
+  real gry = dVdY(ix, iy);
+  if(grx == Fillval() || gry == Fillval()) return Fillval();
+  return michlib::Hypot(grx, gry);
+ }
+
+ real Grad(size_t i) const { return Grad(i % Nx(), i / Nx()); }
 };
 
 class Rect2DData: public BaseData
 {
  std::vector<real> lon, lat;
 
+ static real D(real lon1, real lat1, real lon2, real lat2) { return michlib::GCD(lon1 * M_PI / 180.0, lat1 * M_PI / 180.0, lon2 * M_PI / 180.0, lat2 * M_PI / 180.0) * 6371.0; }
+
  public:
  using BaseData::IsFill;
  using BaseData::V;
@@ -92,4 +151,56 @@ class Rect2DData: public BaseData
 
  real Ix2Lon(size_t ix) const { return lon[ix]; }
  real Iy2Lat(size_t iy) const { return lat[iy]; }
+
+ real dVdX(size_t ix, size_t iy) const
+ {
+  if(IsFill(ix, iy)) return Fillval();
+  if(ix == 0 || IsFill(ix - 1, iy))
+  {
+   if(IsFill(ix + 1, iy))
+    return Fillval();
+   else
+    return (V(ix + 1, iy) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy));
+  }
+  if(ix == Nx() - 1 || IsFill(ix + 1, iy))
+  {
+   if(IsFill(ix - 1, iy))
+    return Fillval();
+   else
+    return (V(ix, iy) - V(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy));
+  }
+  return 0.5 * ((V(ix + 1, iy) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
+                (V(ix, iy) - V(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
+ }
+ real dVdY(size_t ix, size_t iy) const
+ {
+  if(IsFill(ix, iy)) return Fillval();
+
+  if(iy == 0 || IsFill(ix, iy - 1))
+  {
+   if(IsFill(ix, iy + 1))
+    return Fillval();
+   else
+    return (V(ix, iy + 1) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1));
+  }
+  if(iy == Ny() - 1 || IsFill(ix, iy + 1))
+  {
+   if(IsFill(ix, iy - 1))
+    return Fillval();
+   else
+    return (V(ix, iy) - V(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy));
+  }
+  return 0.5 * ((V(ix, iy + 1) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
+                (V(ix, iy) - V(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
+ }
+
+ real Grad(size_t ix, size_t iy) const
+ {
+  real grx = dVdX(ix, iy);
+  real gry = dVdY(ix, iy);
+  if(grx == Fillval() || gry == Fillval()) return Fillval();
+  return michlib::Hypot(grx, gry);
+ }
+
+ real Grad(size_t i) const { return Grad(i % Nx(), i / Nx()); }
 };