NCFuncs::CheckVar now returns a VarPresence enum instead of a bool.

The calculation of derived variables has been moved to a separate function in NCFuncs.
Added calculation of geostrophic velocities based on sea surface height.
The uv and genintfile actions can now use geostrophic velocities.

actions/actiongenintfile.h

@@ -23,6 +23,8 @@ template<class D> MString ActionGenIntFile::DoAction(const CLArgs& args, D& ds)
 
  if(!args.contains("out")) return "Output file name not specified";
 
+ bool mode = args.contains("geostrophic");
+
  auto [tindexes, err] = GetTIndexes(ds, args, pars);
  if(err.Exist()) return err;
  if(tindexes.size() < 10) return "Too few time moments";
@@ -56,7 +58,7 @@ template<class D> MString ActionGenIntFile::DoAction(const CLArgs& args, D& ds)
 
  for(size_t it = 0; it < tindexes.size(); it++)
  {
-  auto data = ReadUV(ds, p, tindexes[it]);
+  auto data = ReadUV(ds, p, tindexes[it], mode);
   if(!data) return "Can't read data";
   if(it == 0)
   {
@@ -91,6 +93,7 @@ template<class D> MString ActionGenIntFile::DoAction(const CLArgs& args, D& ds)
    fw.SetParameter("EndDate", ds.Time(tindexes.back()).ToTString());
    fw.SetParameter("Timestep", ds.Time(tindexes[1]) - ds.Time(tindexes[0]));
    fw.SetParameter("DataID", michlib::UniqueId());
+   fw.SetParameter("Creation command", args.at("_cmdline"));
   }
 
   for(size_t ilat = 0; ilat < data.Ny(); ilat++)

actions/actiontsc.h

@@ -47,7 +47,7 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
   for(const auto& vname: vnames)
   {
    if(used.contains(vname)) return "Duplicate variable " + vname + " in list " + varstring;
-   if(!ds.CheckVar(vname)) return "Variable " + vname + " not exists in this dataset";
+   if(ds.CheckVar(vname) == VarPresence::NONE) return "Variable " + vname + " not exists in this dataset";
    used.insert(vname);
   }
  }

actions/actionuv.h

@@ -106,6 +106,8 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
  pars.UsePrefix("");
  pars.SetParameter("source", args.at("source"));
 
+ bool mode = args.contains("geostrophic");
+
  auto [tindexes, err] = GetTIndexes(ds, args, pars);
  if(err.Exist()) return err;
 
@@ -127,7 +129,7 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
  }
  auto p = sourcepars.get();
 
- auto data = ReadUV(ds, p, tindexes);
+ auto data = ReadUV(ds, p, tindexes, mode);
  if(!data) return "Can't read data";
 
  // Main file

include/actiondep.h

@@ -12,7 +12,7 @@ using michlib::MDateTime;
 #define ADD_ACTION(actclass, actname, suptest, ...) ADD ACTION CLASS: actclass
 #else
 #define ADD_ACTION(actclass, actname, suptest, ...)                                                                                                                                \
- class Action##actclass __VA_OPT__( : protected) __VA_ARGS__                                                                                                                          \
+ class Action##actclass __VA_OPT__( : protected) __VA_ARGS__                                                                                                                       \
  {                                                                                                                                                                                 \
   public:                                                                                                                                                                          \
   static constexpr const char*                 name        = #actname;                                                                                                             \
@@ -151,29 +151,33 @@ template<class D> std::vector<ReadType<D>> Read(const D& data, const std::vector
  return out;
 }
 
-template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, size_t ind)
+template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, size_t ind, bool usegeo)
 {
  using RT = ReadType<D>;
  using UV = UVData<RT>;
  michlib::message("Time: " + data.Time(ind).ToTString());
  std::map<MString, RT> cache;
- bool                  res = false;
+ bool                  res   = false;
+ const MString         uname = usegeo ? "ugeo" : "u";
+ const MString         vname = usegeo ? "vgeo" : "v";
  if constexpr(ReadPSupported<D>)
-  res = data.Read("u", cache, p, ind) && data.Read("v", cache, p, ind);
+  res = data.Read(uname, cache, p, ind) && data.Read(vname, cache, p, ind);
  else if constexpr(ReadSupported<D>)
-  res = data.Read("u", cache, ind) && data.Read("v", cache, ind);
+  res = data.Read(uname, cache, ind) && data.Read(vname, cache, ind);
  if(!res) return UV();
 
- return UV(cache.at("u"), cache.at("v"));
+ return UV(cache.at(uname), cache.at(vname));
 }
 
-template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex)
+template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex, bool usegeo)
 {
- using RT = ReadType<D>;
- using UV = UVData<RT>;
+ using RT            = ReadType<D>;
+ using UV            = UVData<RT>;
+ const MString uname = usegeo ? "ugeo" : "u";
+ const MString vname = usegeo ? "vgeo" : "v";
 
  if(tindex.size() == 1)
-  return ReadUV(data, p, tindex[0]);
+  return ReadUV(data, p, tindex[0], usegeo);
  else
  {
   Averager<UV> out;
@@ -187,12 +191,12 @@ template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters
    std::map<MString, RT> cache;
    bool                  res = false;
    if constexpr(ReadPSupported<D>)
-    res = data.Read("u", cache, p, ind) && data.Read("v", cache, p, ind);
+    res = data.Read(uname, cache, p, ind) && data.Read(vname, cache, p, ind);
    else if constexpr(ReadSupported<D>)
-    res = data.Read("u", cache, ind) && data.Read("v", cache, ind);
+    res = data.Read(uname, cache, ind) && data.Read(vname, cache, ind);
    if(!res) return UV();
 
-   UV dat(cache.at("u"), cache.at("v"));
+   UV dat(cache.at(uname), cache.at(vname));
    if(dat)
     out.Add(std::move(dat));
    else

include/basedata.h

@@ -15,6 +15,13 @@ struct Region
  real lonb, lone, latb, late;
 };
 
+enum class VarPresence
+{
+ NONE,
+ INTERNAL,
+ DERIVED
+};
+
 class BaseData
 {
  protected:

include/layereddata.h

@@ -143,6 +143,8 @@ class LayeredData: public NCFuncs
 
  real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; }
 
+ real Depth(const BaseParameters* ip) const { return Depth(dynamic_cast<const struct Parameters*>(ip)->layer); }
+
  real Lon(size_t ix) const { return isOk() ? (lonb + ix * lonstep) : -1000.0; }
 
  real Lat(size_t iy) const { return isOk() ? (latb + iy * latstep) : -1000.0; }
@@ -172,7 +174,7 @@ class LayeredData: public NCFuncs
   // clang-format on
  }
 
- bool CheckVar(const MString& vname) const
+ VarPresence CheckVar(const MString& vname) const
  {
   return NCFuncs::CheckVar(vname, [this](const MString& vn) { return HaveVar(vn); });
  }

include/ncfuncs.h

@@ -1,6 +1,8 @@
 #pragma once
 #include "DataAdapters/ncfilealt.h"
+#include "basedata.h"
 #include "mdatetime.h"
+#include <map>
 #include <set>
 #include <tuple>
 
@@ -25,17 +27,197 @@ class NCFuncs
  static CoordNames                          GetCNames(const NCFileA& nc);
  static CoordNames                          GetDNames(const NCFileA& nc);
  static bool                                HaveVar(const NCFileA& nc, const MString& vname);
- template<class HV> static bool             CheckVar(const MString& vname, HV hv)
+ template<class HV> static VarPresence      CheckVar(const MString& vname, HV hv)
  {
   if(!hv(vname))
   {
-   bool varexist = false;
-   if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = true;
-   if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = true;
-   if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = true;
-   if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = true;
-   if(!varexist) return false;
+   auto varexist = VarPresence::NONE;
+   if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = VarPresence::DERIVED;
+   if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = VarPresence::DERIVED;
+   if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = VarPresence::DERIVED;
+   if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = VarPresence::DERIVED;
+   if((vname == "ugeo" || vname == "vgeo") && hv("ssh")) varexist = VarPresence::DERIVED;
+   return varexist;
   }
-  return true;
+  return VarPresence::INTERNAL;
+ }
+
+ template<class D> static bool TransformationRead(const D* data, const MString& vname, std::map<MString, typename D::Data>& cache, const BaseParameters* ip, size_t i)
+ {
+  if constexpr(requires(const D* d, const BaseParameters* p) {
+                {
+                 d->Depth(p)
+                } -> std::convertible_to<real>;
+               })
+  {
+   real depth = data->Depth(ip);
+   // ptemp from temp and sal
+   if(vname == "ptemp")
+   {
+    if(!(data->Read("temp", cache, ip, i) && data->Read("sal", cache, ip, i))) return false;
+    cache[vname]      = cache.at("temp").CopyGrid();
+    auto&       ptemp = cache.at(vname);
+    const auto& temp  = cache.at("temp");
+    const auto& sal   = cache.at("sal");
+    ptemp.SetUnit("degrees_C");
+    for(size_t ind = 0; ind < ptemp.N(); ind++)
+    {
+     if(temp.IsFill(ind) || sal.IsFill(ind))
+      ptemp.V(ind) = ptemp.Fillval();
+     else
+      ptemp.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), depth, ptemp.Lon(ind), ptemp.Lat(ind));
+    }
+    return true;
+   }
+   // temp from ptemp and sal
+   if(vname == "temp")
+   {
+    if(!(data->Read("ptemp", cache, ip, i) && data->Read("sal", cache, ip, i))) return false;
+    cache[vname]      = cache.at("ptemp").CopyGrid();
+    auto&       temp  = cache.at(vname);
+    const auto& ptemp = cache.at("ptemp");
+    const auto& sal   = cache.at("sal");
+    temp.SetUnit("degrees_C");
+    for(size_t ind = 0; ind < temp.N(); ind++)
+    {
+     if(ptemp.IsFill(ind) || sal.IsFill(ind))
+      temp.V(ind) = temp.Fillval();
+     else
+      temp.V(ind) = PTemp2Temp(ptemp.V(ind), sal.V(ind), depth, temp.Lon(ind), temp.Lat(ind));
+    }
+    return true;
+   }
+   // pdens from temp and sal
+   if(vname == "pdens")
+   {
+    bool tempispot = data->CheckVar("ptemp") == VarPresence::INTERNAL;
+    if(!(data->Read(tempispot ? "ptemp" : "temp", cache, ip, i) && data->Read("sal", cache, ip, i))) return false;
+    cache[vname]      = cache.at("sal").CopyGrid();
+    auto&       pdens = cache.at(vname);
+    const auto& temp  = cache.at(tempispot ? "ptemp" : "temp");
+    const auto& sal   = cache.at("sal");
+    pdens.SetUnit("kg m-3");
+    for(size_t ind = 0; ind < pdens.N(); ind++)
+    {
+     if(temp.IsFill(ind) || sal.IsFill(ind))
+      pdens.V(ind) = pdens.Fillval();
+     else
+      pdens.V(ind) =
+          tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), depth, pdens.Lon(ind), pdens.Lat(ind)) : Temp2PDens(temp.V(ind), sal.V(ind), depth, pdens.Lon(ind), pdens.Lat(ind));
+    }
+    return true;
+   }
+  }
+
+  // U and U2 from u and v
+  if(vname == "U" || vname == "U2")
+  {
+   bool square = vname == "U2";
+   if(!(data->Read("u", cache, ip, i) && data->Read("v", cache, ip, i))) return false;
+   cache[vname]  = cache.at("u").CopyGrid();
+   auto&       U = cache.at(vname);
+   const auto& u = cache.at("u");
+   const auto& v = cache.at("v");
+   if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
+   for(size_t ind = 0; ind < U.N(); ind++)
+   {
+    if(u.IsFill(ind) || v.IsFill(ind))
+     U.V(ind) = U.Fillval();
+    else
+     U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
+   }
+   return true;
+  }
+  // ugeo
+  if(vname == "ugeo")
+  {
+   if(!data->Read("ssh", cache, ip, i)) return false;
+   const auto& ssh = cache.at("ssh");
+   const real  mul = (ssh.Unit() == "cm") ? 1.0 : ((ssh.Unit() == "m") ? 100.0 : -1.0);
+   if(mul < 0.0) return false; // Unknown units
+
+   cache[vname] = cache.at("ssh").CopyGrid();
+   auto& ugeo   = cache.at(vname);
+   ugeo.SetUnit("cm/s");
+
+   for(size_t iy = 0; iy < ssh.Ny(); iy++)
+   {
+    const real phi   = michlib::M_PI * ssh.Iy2Lat(iy) / 180.0;
+    const real sphi  = michlib::Sin(phi);
+    const real s2phi = michlib::Sin(2.0 * phi);
+    const real g     = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi);                    // g in cm/s
+    const real f     = 2.0 * 7.29211 * sphi;                                                                    // Coriolis parameter multiplied by 100'000
+    const real db = (iy == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Iy2Lat(iy) - ssh.Iy2Lat(iy - 1)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter
+    const real du =
+        (iy == ssh.Ny() - 1) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Iy2Lat(iy + 1) - ssh.Iy2Lat(iy)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter
+
+    for(size_t ix = 0; ix < ssh.Nx(); ix++)
+    {
+     const bool bok = iy != 0 && !ssh.IsFill(ix, iy - 1) && !ssh.IsFill(ix, iy);
+     const bool uok = iy != ssh.Ny() - 1 && !ssh.IsFill(ix, iy + 1) && !ssh.IsFill(ix, iy);
+     if(!(bok || uok))
+      ugeo(ix, iy) = ugeo.Fillval();
+     else if(!bok)
+      ugeo(ix, iy) = -g * mul * (ssh(ix, iy + 1) - ssh(ix, iy)) / (f * du);
+     else if(!uok)
+      ugeo(ix, iy) = -g * mul * (ssh(ix, iy) - ssh(ix, iy - 1)) / (f * db);
+     else
+     {
+      const real bder = -g * mul * (ssh(ix, iy) - ssh(ix, iy - 1)) / (f * db);
+      const real uder = -g * mul * (ssh(ix, iy + 1) - ssh(ix, iy)) / (f * du);
+      ugeo(ix, iy)    = 0.5 * (bder + uder);
+     }
+    }
+   }
+   return true;
+  }
+  // vgeo
+  if(vname == "vgeo")
+  {
+   if(!data->Read("ssh", cache, ip, i)) return false;
+   const auto& ssh = cache.at("ssh");
+   const real  mul = (ssh.Unit() == "cm") ? 1.0 : ((ssh.Unit() == "m") ? 100.0 : -1.0);
+   if(mul < 0.0) return false; // Unknown units
+
+   cache[vname] = cache.at("ssh").CopyGrid();
+   auto& vgeo   = cache.at(vname);
+   vgeo.SetUnit("cm/s");
+
+   for(size_t iy = 0; iy < ssh.Ny(); iy++)
+   {
+    const real phi   = michlib::M_PI * ssh.Iy2Lat(iy) / 180.0;
+    const real sphi  = michlib::Sin(phi);
+    const real cphi  = michlib::Cos(phi);
+    const real s2phi = michlib::Sin(2.0 * phi);
+    const real g     = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi); // g in cm/s
+    const real f     = 2.0 * 7.29211 * sphi;                                                 // Coriolis parameter multiplied by 100'000
+
+    for(size_t ix = 0; ix < ssh.Nx(); ix++)
+    {
+     const real dl =
+         (ix == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Ix2Lon(ix) - ssh.Ix2Lon(ix - 1)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter
+     const real dr = (ix == ssh.Nx() - 1)
+                         ? 0.0
+                         : (michlib::M_PI * 6371.0 * (ssh.Ix2Lon(ix + 1) - ssh.Ix2Lon(ix)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter
+
+     const bool lok = ix != 0 && !ssh.IsFill(ix - 1, iy) && !ssh.IsFill(ix, iy);
+     const bool rok = ix != ssh.Nx() - 1 && !ssh.IsFill(ix + 1, iy) && !ssh.IsFill(ix, iy);
+     if(!(lok || rok))
+      vgeo(ix, iy) = vgeo.Fillval();
+     else if(!lok)
+      vgeo(ix, iy) = g * mul * (ssh(ix + 1, iy) - ssh(ix, iy)) / (f * dr);
+     else if(!rok)
+      vgeo(ix, iy) = g * mul * (ssh(ix, iy) - ssh(ix - 1, iy)) / (f * dl);
+     else
+     {
+      const real lder = g * mul * (ssh(ix, iy) - ssh(ix - 1, iy)) / (f * dl);
+      const real rder = g * mul * (ssh(ix + 1, iy) - ssh(ix, iy)) / (f * dr);
+      vgeo(ix, iy)    = 0.5 * (lder + rder);
+     }
+    }
+   }
+   return true;
+  }
+  return false;
  }
 };

include/traits.h

@@ -167,3 +167,33 @@ concept ReadIs1DGeoArray = requires {
   std::declval<ReadType<T>>().Lat(0)
  } -> std::convertible_to<real>;
 };
+
+template<class T>
+concept ReadIs2DXYRectArray = requires {
+ {
+  std::declval<ReadType<T>>().Ix2X(0)
+ } -> std::convertible_to<real>;
+ {
+  std::declval<ReadType<T>>().Iy2Y(0)
+ } -> std::convertible_to<real>;
+};
+
+template<class T>
+concept ReadIs2DXYArray = requires {
+ {
+  std::declval<ReadType<T>>().X(0, 0)
+ } -> std::convertible_to<real>;
+ {
+  std::declval<ReadType<T>>().Y(0, 0)
+ } -> std::convertible_to<real>;
+};
+
+template<class T>
+concept ReadIs1DArray = requires {
+ {
+  std::declval<ReadType<T>>().X(0)
+ } -> std::convertible_to<real>;
+ {
+  std::declval<ReadType<T>>().Y(0)
+ } -> std::convertible_to<real>;
+};

sources/AVISOLOCAL.cpp

@@ -115,28 +115,7 @@ bool AVISOLOCALData::Read(const MString& vname, std::map<MString, AVISOLOCALData
  }
 
  if(!name.Exist()) // Conversion read
- {
-  // U and U2 from u and v
-  if(vname == "U" || vname == "U2")
-  {
-   bool square = vname == "U2";
-   if(!(Read("u", cache, ip, i) && Read("v", cache, ip, i))) return false;
-   cache[vname]  = cache.at("u").CopyGrid();
-   auto&       U = cache.at(vname);
-   const auto& u = cache.at("u");
-   const auto& v = cache.at("v");
-   if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
-   for(size_t ind = 0; ind < U.N(); ind++)
-   {
-    if(u.IsFill(ind) || v.IsFill(ind))
-     U.V(ind) = U.Fillval();
-    else
-     U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
-   }
-   return true;
-  }
-  return false;
- }
+  return TransformationRead(this, vname, cache, ip, i);
 
  // Direct read
  Data data;

sources/AVISOLOCAL.h

@@ -53,7 +53,7 @@ class AVISOLOCALData: public NCFuncs
 
  explicit operator bool() const { return times.size() > 0; }
 
- bool CheckVar(const MString& vname) const
+ VarPresence CheckVar(const MString& vname) const
  {
   NCFileA nc;
   nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc");

sources/BINFILE.h

@@ -42,7 +42,12 @@ class BINFILEData
 
  explicit operator bool() const { return times.size() > 0; }
 
- bool CheckVar(const MString& vname) const { return vname == "u" || vname == "v" || vname == "U" || vname == "U2"; }
+ VarPresence CheckVar(const MString& vname) const
+ {
+  if(vname == "u" || vname == "v") return VarPresence::INTERNAL;
+  if(vname == "U" || vname == "U2") return VarPresence::DERIVED;
+  return VarPresence::NONE;
+ }
 
  bool Read(const MString& vname, std::map<MString, Data>& cache, size_t i) const;
 };

sources/MODISBINLOCAL.h

@@ -54,11 +54,12 @@ class MODISBINLOCALData
 
  std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const;
 
- bool CheckVar(const MString& vname) const
+ VarPresence CheckVar(const MString& vname) const
  {
-  if(dataset.empty()) return false;
+  if(dataset.empty()) return VarPresence::NONE;
   bool ischl = dataset.front() == "A_CHL" || dataset.front() == "T_CHL";
-  return (vname == "chl" && ischl) || (vname == "temp" && !ischl);
+  if((vname == "chl" && ischl) || (vname == "temp" && !ischl)) return VarPresence::INTERNAL;
+  return VarPresence::NONE;
  }
 
  bool isOk() const { return times.size() > 0; }

sources/VYLET.h

@@ -71,14 +71,14 @@ class VYLETData
   return vars;
  }
 
- bool CheckVar(const MString& vname) const
+ VarPresence CheckVar(const MString& vname) const
  {
   std::set<MString> vars{"vylD", "vylNx", "vylNy", "vylRx", "vylRy", "vylEx", "vylEy", "vylPhip", "vylPhim", "vylPhit", "vylAngle", "vylLen"};
   if(HasLyap()) vars.insert("vylL");
   if(HasLyap()) vars.insert("vyldS");
   if(HasTime()) vars.insert("vylT");
   if(HasTime()) vars.insert("vylTmask");
-  return vars.contains(vname);
+  return vars.contains(vname) ? VarPresence::DERIVED : VarPresence::NONE;
  }
 
  bool Read(const MString& vname, std::map<MString, Data>& cache, size_t tind) const;

src/layereddata.cpp

@@ -186,84 +186,7 @@ bool LayeredData::Read(const MString& vname, std::map<MString, LayeredData::Data
  auto p               = dynamic_cast<const struct Parameters*>(ip);
  auto [name, id, tid] = VarNameLoc(vname, times[i]);
  if(!name.Exist()) // Conversion read
- {
-  // ptemp from temp and sal
-  if(vname == "ptemp")
-  {
-   if(!(Read("temp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
-   cache[vname]      = cache.at("temp").CopyGrid();
-   auto&       ptemp = cache.at(vname);
-   const auto& temp  = cache.at("temp");
-   const auto& sal   = cache.at("sal");
-   ptemp.SetUnit("degrees_C");
-   for(size_t ind = 0; ind < ptemp.N(); ind++)
-   {
-    if(temp.IsFill(ind) || sal.IsFill(ind))
-     ptemp.V(ind) = ptemp.Fillval();
-    else
-     ptemp.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), Depth(p->layer), ptemp.Lon(ind), ptemp.Lat(ind));
-   }
-   return true;
-  }
-  // temp from ptemp and sal
-  if(vname == "temp")
-  {
-   if(!(Read("ptemp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
-   cache[vname]      = cache.at("ptemp").CopyGrid();
-   auto&       temp  = cache.at(vname);
-   const auto& ptemp = cache.at("ptemp");
-   const auto& sal   = cache.at("sal");
-   temp.SetUnit("degrees_C");
-   for(size_t ind = 0; ind < temp.N(); ind++)
-   {
-    if(ptemp.IsFill(ind) || sal.IsFill(ind))
-     temp.V(ind) = temp.Fillval();
-    else
-     temp.V(ind) = PTemp2Temp(ptemp.V(ind), sal.V(ind), Depth(p->layer), temp.Lon(ind), temp.Lat(ind));
-   }
-   return true;
-  }
-  // pdens from temp and sal
-  if(vname == "pdens")
-  {
-   bool tempispot = HaveVar("ptemp");
-   if(!(Read(tempispot ? "ptemp" : "temp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
-   cache[vname]      = cache.at("sal").CopyGrid();
-   auto&       pdens = cache.at(vname);
-   const auto& temp  = cache.at(tempispot ? "ptemp" : "temp");
-   const auto& sal   = cache.at("sal");
-   pdens.SetUnit("kg m-3");
-   for(size_t ind = 0; ind < pdens.N(); ind++)
-   {
-    if(temp.IsFill(ind) || sal.IsFill(ind))
-     pdens.V(ind) = pdens.Fillval();
-    else
-     pdens.V(ind) = tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), pdens.Lon(ind), pdens.Lat(ind))
-                              : Temp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), pdens.Lon(ind), pdens.Lat(ind));
-   }
-   return true;
-  }
-  // U and U2 from u and v
-  if(vname == "U" || vname == "U2")
-  {
-   bool square = vname == "U2";
-   if(!(Read("u", cache, ip, i) && Read("v", cache, ip, i))) return false;
-   cache[vname]  = cache.at("u").CopyGrid();
-   auto&       U = cache.at(vname);
-   const auto& u = cache.at("u");
-   const auto& v = cache.at("v");
-   if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
-   for(size_t ind = 0; ind < U.N(); ind++)
-   {
-    if(u.IsFill(ind) || v.IsFill(ind))
-     U.V(ind) = U.Fillval();
-    else
-     U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
-   }
-   return true;
-  }
-  return false;
- }
+ return TransformationRead(this, vname, cache, ip, i);
 
  // Direct read
  bool nodepth = false;

src/ncfuncs.cpp

@@ -79,6 +79,8 @@ void NCFuncs::GetVars(const NCFileA& nc, std::set<MString>& vars)
 
  if(vars.contains("u") && vars.contains("v")) vars.emplace("U");
  if(vars.contains("u") && vars.contains("v")) vars.emplace("U2");
+ if(vars.contains("ssh")) vars.emplace("ugeo");
+ if(vars.contains("ssh")) vars.emplace("vgeo");
 }
 
 std::tuple<MDateTime, time_t, bool> NCFuncs::Refdate(const MString& refdate)
@@ -155,6 +157,8 @@ MString NCFuncs::Name2LongName(const MString& name)
  if(name == "u") return "X-velocity";
  if(name == "v") return "Y-velocity";
  if(name == "w") return "Z-velocity";
+ if(name == "ugeo") return "Geostrophic eastward velocity from ssh";
+ if(name == "vgeo") return "Geostrophic northward velocity from ssh";
  if(name == "chl") return "Concentration of chlorophyll";
  if(name == "NO3") return "Concentration of nitrates";
  if(name == "prprod") return "Primary production";