From bf659b1b81156809c982e93c27ca5301de52279e Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Fri, 29 Sep 2023 16:02:55 +1000 Subject: [PATCH] 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 | 5 +- actions/actiontsc.h | 2 +- actions/actionuv.h | 4 +- include/actiondep.h | 30 +++--- include/basedata.h | 7 ++ include/layereddata.h | 4 +- include/ncfuncs.h | 198 +++++++++++++++++++++++++++++++++++-- sources/AVISOLOCAL.cpp | 23 +---- sources/AVISOLOCAL.h | 2 +- sources/BINFILE.h | 7 +- sources/MODISBINLOCAL.h | 7 +- sources/VYLET.h | 4 +- src/layereddata.cpp | 79 +-------------- src/ncfuncs.cpp | 4 + 14 files changed, 244 insertions(+), 132 deletions(-) diff --git a/actions/actiongenintfile.h b/actions/actiongenintfile.h index 7b04a21..42a09d7 100644 --- a/actions/actiongenintfile.h +++ b/actions/actiongenintfile.h @@ -23,6 +23,8 @@ template 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 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 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++) diff --git a/actions/actiontsc.h b/actions/actiontsc.h index a780c39..fabda5a 100644 --- a/actions/actiontsc.h +++ b/actions/actiontsc.h @@ -47,7 +47,7 @@ template 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); } } diff --git a/actions/actionuv.h b/actions/actionuv.h index 01c96f7..62dbb79 100644 --- a/actions/actionuv.h +++ b/actions/actionuv.h @@ -106,6 +106,8 @@ template 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 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 diff --git a/include/actiondep.h b/include/actiondep.h index b1f7e14..2cc8367 100644 --- a/include/actiondep.h +++ b/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 std::vector> Read(const D& data, const std::vector return out; } -template UVData> ReadUV(const D& data, const BaseParameters* p, size_t ind) +template UVData> ReadUV(const D& data, const BaseParameters* p, size_t ind, bool usegeo) { using RT = ReadType; using UV = UVData; michlib::message("Time: " + data.Time(ind).ToTString()); std::map cache; - bool res = false; + bool res = false; + const MString uname = usegeo ? "ugeo" : "u"; + const MString vname = usegeo ? "vgeo" : "v"; if constexpr(ReadPSupported) - 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) - 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 UVData> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex) +template UVData> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex, bool usegeo) { - using RT = ReadType; - using UV = UVData; + using RT = ReadType; + using UV = UVData; + 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 out; @@ -187,12 +191,12 @@ template UVData> ReadUV(const D& data, const BaseParameters std::map cache; bool res = false; if constexpr(ReadPSupported) - 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) - 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 diff --git a/include/basedata.h b/include/basedata.h index a9a03e1..5b8948b 100644 --- a/include/basedata.h +++ b/include/basedata.h @@ -15,6 +15,13 @@ struct Region real lonb, lone, latb, late; }; +enum class VarPresence +{ + NONE, + INTERNAL, + DERIVED +}; + class BaseData { protected: diff --git a/include/layereddata.h b/include/layereddata.h index 9e5b984..b97762c 100644 --- a/include/layereddata.h +++ b/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(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); }); } diff --git a/include/ncfuncs.h b/include/ncfuncs.h index 2dcf444..207f4af 100644 --- a/include/ncfuncs.h +++ b/include/ncfuncs.h @@ -1,6 +1,8 @@ #pragma once #include "DataAdapters/ncfilealt.h" +#include "basedata.h" #include "mdatetime.h" +#include #include #include @@ -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 static bool CheckVar(const MString& vname, HV hv) + template 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 static bool TransformationRead(const D* data, const MString& vname, std::map& cache, const BaseParameters* ip, size_t i) + { + if constexpr(requires(const D* d, const BaseParameters* p) { + { + d->Depth(p) + } -> std::convertible_to; + }) + { + 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; } }; diff --git a/sources/AVISOLOCAL.cpp b/sources/AVISOLOCAL.cpp index 1847d58..a61fd1c 100644 --- a/sources/AVISOLOCAL.cpp +++ b/sources/AVISOLOCAL.cpp @@ -115,28 +115,7 @@ bool AVISOLOCALData::Read(const MString& vname, std::map 0; } - bool CheckVar(const MString& vname) const + VarPresence CheckVar(const MString& vname) const { NCFileA nc; nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc"); diff --git a/sources/BINFILE.h b/sources/BINFILE.h index ad9d7b9..eb3cc18 100644 --- a/sources/BINFILE.h +++ b/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& cache, size_t i) const; }; diff --git a/sources/MODISBINLOCAL.h b/sources/MODISBINLOCAL.h index e5ed801..ca2ab07 100644 --- a/sources/MODISBINLOCAL.h +++ b/sources/MODISBINLOCAL.h @@ -54,11 +54,12 @@ class MODISBINLOCALData std::pair 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; } diff --git a/sources/VYLET.h b/sources/VYLET.h index 9c9089a..ffe7b90 100644 --- a/sources/VYLET.h +++ b/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 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& cache, size_t tind) const; diff --git a/src/layereddata.cpp b/src/layereddata.cpp index f15ab7b..dc9f8f6 100644 --- a/src/layereddata.cpp +++ b/src/layereddata.cpp @@ -186,84 +186,7 @@ bool LayeredData::Read(const MString& vname, std::map(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; diff --git a/src/ncfuncs.cpp b/src/ncfuncs.cpp index 3bb2745..5ab0d3b 100644 --- a/src/ncfuncs.cpp +++ b/src/ncfuncs.cpp @@ -79,6 +79,8 @@ void NCFuncs::GetVars(const NCFileA& nc, std::set& 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 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";