From 4f6782ea098dd9cad35df7eaa8a720f9316a9621 Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Tue, 7 Mar 2023 15:03:28 +1000 Subject: [PATCH] More variables. Update NEMO code. --- include/NEMO.h | 199 ++++++++++++++++++++++++++-------------------- include/vartype.h | 33 +++++--- src/actiontsc.cpp | 2 +- src/odm.cpp | 6 +- 4 files changed, 139 insertions(+), 101 deletions(-) diff --git a/include/NEMO.h b/include/NEMO.h index 273ed26..614bd5d 100644 --- a/include/NEMO.h +++ b/include/NEMO.h @@ -27,7 +27,7 @@ class NEMOData TYPE_NRT6 }; - NCFileA nc, nct, ncs; + std::vector nc; size_t xb, yb, xe, ye, layer; std::vector depths; std::vector times; @@ -73,15 +73,22 @@ class NEMOData public: using Data = Simple2DData; - Data ReadVar(const NCFileA& f, const MString& name, size_t i) const + private: + template Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth) const { - using DataType = int2; - constexpr DataType fill = -32767; - real unitmul = 1.0; + real unitmul = 1.0; + DataType fill; + real offset = 0.0, scale = 1.0; - auto offset = f.A(name, "add_offset"); - auto scale = f.A(name, "scale_factor"); - if(!offset || !scale) return Data(); + { + auto a_fill = f.A(name, "_FillValue"); + auto a_offset = f.A(name, "add_offset"); + auto a_scale = f.A(name, "scale_factor"); + if(!a_fill) return Data(); + fill = a_fill; + if(a_offset) offset = a_offset; + if(a_scale) scale = a_scale; + } auto unit = f.A(name, "units"); if(unit && unit.Get() == "m s-1") unitmul = 100.0; @@ -90,7 +97,8 @@ class NEMOData if(xb < xe) { - auto var = f.V(name, {"longitude", xb, xe - xb + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); + auto var = nodepth ? f.V(name, {"longitude", xb, xe - xb + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}) + : f.V(name, {"longitude", xb, xe - xb + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); if(!var) return Data(); if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data(); @@ -103,8 +111,10 @@ class NEMOData } else { - auto var1 = f.V(name, {"longitude", xb}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); - auto var2 = f.V(name, {"longitude", 0, xe + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); + auto var1 = nodepth ? f.V(name, {"longitude", xb}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}) + : f.V(name, {"longitude", xb}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); + auto var2 = nodepth ? f.V(name, {"longitude", 0, xe + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}) + : f.V(name, {"longitude", 0, xe + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); if(!(var1 && var2)) return Data(); if((var1.DimLen(0) + var2.DimLen(0)) != data.Nx() || var1.DimLen(1) != data.Ny() || var2.DimLen(1) != data.Ny()) return Data(); for(size_t ix = 0; ix < var1.DimLen(0); ix++) @@ -123,15 +133,35 @@ class NEMOData return data; } - Data ReadVar(const NCFileA& f, const MString& name, const std::vector& tindex) const + public: + Data ReadVar(const MString& name, size_t i) const + { + if(!isOk()) return Data(); + bool nodepth = false; + + for(const auto& f: nc) + { + auto head = f.Header(); + for(const auto& v: head.Variables()) + if(v.Name() == name) + { + if(v.Dimensions().size() == 3) nodepth = true; + if(v.Type().Id() == NC_SHORT) return ReadVarRaw(f, name, i, nodepth); + if(v.Type().Id() == NC_FLOAT) return ReadVarRaw(f, name, i, nodepth); + } + } + return Data(); + } + + Data ReadVar(const MString& name, const std::vector& tindex) const { Data out; - if(tindex.size() == 0) return out; + if(tindex.size() == 0 || !isOk()) return out; std::vector count; for(size_t i = 0; i < tindex.size(); i++) { - Data dat = ReadVar(f, name, tindex[i]); + Data dat = ReadVar(name, tindex[i]); if(!dat) return Data(); out.Add(dat, count); } @@ -145,76 +175,69 @@ class NEMOData // TODO: RetVal bool Open(const MString& stype, const MString& cred, const MString& proxyurl = "") { + nc.clear(); if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl); - MString url, urlt, urls; if(stype == "DT") { - url = "https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m"; + NCFileA newnc; + newnc.Reset("https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m"); + if(!newnc) return false; + nc.push_back(std::move(newnc)); type = TYPE_DT; } if(stype == "NRT") { - url = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024"; - type = TYPE_NRT; - } - if(stype == "NRT6") - { - url = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-uovo"; - urlt = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-thetao"; - urls = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-so"; - type = TYPE_NRT6; - } - nc.Reset(url); - if(!nc) return false; - - { - auto nx_ = nc.D("longitude"); - if(!nx_) - { - nc.Reset(); - return false; - } - nx = nx_; - } - - if(urlt.Exist()) - { - nct.Reset(urlt); - if(!nct) + std::vector urls{"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-cur_anfc_0.083deg_P1D-m", + "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-thetao_anfc_0.083deg_P1D-m", + "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-so_anfc_0.083deg_P1D-m", + "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_anfc_0.083deg_P1D-m", + "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-wcur_anfc_0.083deg_P1D-m"}; + for(const auto& url: urls) { - nc.Reset(); - return false; + NCFileA newnc; + newnc.Reset(url); + if(!newnc) + { + nc.clear(); + return false; + } + nc.push_back(std::move(newnc)); } + type = TYPE_NRT; } - if(urls.Exist()) + if(stype == "NRT6") { - ncs.Reset(urls); - if(!ncs) + std::vector urls{"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-cur_anfc_0.083deg_PT6H-i", + "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-thetao_anfc_0.083deg_PT6H-i", + "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-so_anfc_0.083deg_PT6H-i"}; + for(const auto& url: urls) { - nc.Reset(); - nct.Reset(); - return false; + NCFileA newnc; + newnc.Reset(url); + if(!newnc) + { + nc.clear(); + return false; + } + nc.push_back(std::move(newnc)); } + type = TYPE_NRT6; } - auto rdepths = nc.V("depth"); + auto rdepths = nc[0].V("depth"); if(!rdepths) { - nc.Reset(); - nct.Reset(); - ncs.Reset(); + nc.clear(); return false; } depths.resize(rdepths.DimLen(0)); for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i); - auto timeD = nc.V("time"); - auto timeF = nc.V("time"); + auto timeD = nc[0].V("time"); + auto timeF = nc[0].V("time"); if(!(timeD || timeF)) { - nc.Reset(); - nct.Reset(); - ncs.Reset(); + nc.clear(); return false; } MDateTime refdate("1950-01-01"); @@ -253,7 +276,7 @@ class NEMOData bool SetRegion(real lonbin, real lonein, real latb, real late, real depth) { - if(!nc) return false; + if(!isOk()) return false; real lonb = ToGeoDomainNeg(lonbin), lone = ToGeoDomainNeg(lonein); size_t xb_, xe_, yb_, ye_, layer_ = 0; @@ -286,28 +309,28 @@ class NEMOData return true; } - real Lonb() const { return nc ? (-180.0 + xb / 12.0) : -1000.0; } - real Lone() const { return nc ? (-180.0 + xe / 12.0) : -1000.0; } - real Latb() const { return nc ? (-80.0 + yb / 12.0) : -1000.0; } - real Late() const { return nc ? (-80.0 + ye / 12.0) : -1000.0; } - size_t Xb() const { return nc ? xb : 0; } - size_t Xe() const { return nc ? xe : 0; } - size_t Yb() const { return nc ? yb : 0; } - size_t Ye() const { return nc ? ye : 0; } - size_t Nx() const { return nc ? ((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1)) : 0; } - size_t Ny() const { return nc ? (ye - yb + 1) : 0; } - real Depth(size_t l) const { return nc ? depths[l] : -1000.0; } + real Lonb() const { return isOk() ? (-180.0 + xb / 12.0) : -1000.0; } + real Lone() const { return isOk() ? (-180.0 + xe / 12.0) : -1000.0; } + real Latb() const { return isOk() ? (-80.0 + yb / 12.0) : -1000.0; } + real Late() const { return isOk() ? (-80.0 + ye / 12.0) : -1000.0; } + size_t Xb() const { return isOk() ? xb : 0; } + size_t Xe() const { return isOk() ? xe : 0; } + size_t Yb() const { return isOk() ? yb : 0; } + size_t Ye() const { return isOk() ? ye : 0; } + size_t Nx() const { return isOk() ? ((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1)) : 0; } + size_t Ny() const { return isOk() ? (ye - yb + 1) : 0; } + real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; } real Depth() const { return Depth(layer); } size_t Layer() const { return layer; } - time_t Timestep() const { return nc ? (times[1] - times[0]) : 0; } + time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; } MDateTime Time(size_t i) const { - if((!nc) || i >= times.size()) return MDateTime(); + if(!isOk() || i >= times.size()) return MDateTime(); return times[i]; } - bool isOk() const { return nc; } - explicit operator bool() const { return nc; } + bool isOk() const { return nc.size() > 0; } + explicit operator bool() const { return nc.size() > 0; } size_t NDepths() const { return depths.size(); } size_t NTimes() const { return times.size(); } @@ -330,14 +353,20 @@ class NEMOData static real Fillval() { return Data::Fillval(); } }; -template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar(nc, "uo", it); } -template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar(nc, "vo", it); } -template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar(nct ? nct : nc, "thetao", it); } -template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar(ncs ? ncs : nc, "so", it); } - -template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar(nc, "uo", tindex); } -template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar(nc, "vo", tindex); } -template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar(nct ? nct : nc, "thetao", tindex); } -template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar(ncs ? ncs : nc, "so", tindex); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar("uo", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar("vo", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar("thetao", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar("so", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar("mlotst", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar("zos", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar("wo", it); } + +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar("uo", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar("vo", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar("thetao", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar("so", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar("mlotst", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar("zos", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar("wo", tindex); } #endif diff --git a/include/vartype.h b/include/vartype.h index 6232973..dc2ecc5 100644 --- a/include/vartype.h +++ b/include/vartype.h @@ -15,7 +15,10 @@ enum class Vartype V, TEMP, SAL, - CHL + CHL, + MLD, + SSH, + W }; template struct VartypeCarrier @@ -23,19 +26,19 @@ template struct VartypeCarrier static const Vartype vt = vt_; }; -using VartypeUnion = std::variant, VartypeCarrier, VartypeCarrier, VartypeCarrier, - VartypeCarrier, VartypeCarrier>; +using VartypeUnion = + std::variant, VartypeCarrier, VartypeCarrier, VartypeCarrier, VartypeCarrier, + VartypeCarrier, VartypeCarrier, VartypeCarrier, VartypeCarrier>; template -concept HasVar = requires(T t) -{ - { - t.template Read(0)(0, 0) - } -> std::convertible_to; - { - t.template Read(std::vector())(0, 0) - } -> std::convertible_to; -}; +concept HasVar = requires(T t) { + { + t.template Read(0)(0, 0) + } -> std::convertible_to; + { + t.template Read(std::vector())(0, 0) + } -> std::convertible_to; + }; } // namespace vartype template static constexpr bool isDataSupported = vartype::HasVar; @@ -59,6 +62,9 @@ class VarType: public vartype::VartypeUnion if(str == "t" || str == "temp" || str == "temperature") *this = vartype::VartypeCarrier(); if(str == "s" || str == "sal" || str == "salinity") *this = vartype::VartypeCarrier(); if(str == "c" || str == "chl" || str == "chlorophyll") *this = vartype::VartypeCarrier(); + if(str == "mld") *this = vartype::VartypeCarrier(); + if(str == "ssh") *this = vartype::VartypeCarrier(); + if(str == "w") *this = vartype::VartypeCarrier(); } MString Name() const @@ -71,6 +77,9 @@ class VarType: public vartype::VartypeUnion case(vartype::Vartype::TEMP): return "Temperature"; case(vartype::Vartype::SAL): return "Salinity"; case(vartype::Vartype::CHL): return "Chlorophyll"; + case(vartype::Vartype::MLD): return "Mixed layer depth"; + case(vartype::Vartype::SSH): return "Sea surface height"; + case(vartype::Vartype::W): return "W"; } return "none"; } diff --git a/src/actiontsc.cpp b/src/actiontsc.cpp index 2ca7364..40bd79e 100644 --- a/src/actiontsc.cpp +++ b/src/actiontsc.cpp @@ -89,7 +89,7 @@ int actiontsc(const CLArgs& args) } BFileW fw; - MString name = args.at("out"); + MString name = args.contains("out") ? args.at("out") : ""; if(!name.Exist()) name = "out.bin"; fw.Create(name, 3); diff --git a/src/odm.cpp b/src/odm.cpp index d6e31ec..cba317f 100644 --- a/src/odm.cpp +++ b/src/odm.cpp @@ -11,14 +11,14 @@ inline void Usage(const MString& arg0) message(" dataset. Can be DT, NRT or NRT6. Default: DT"); message(" Keys for action=tsc. Get temperature, salinity, chlorofill from dataset."); message(" source. Required. May be: NEMO"); - message(" var. Required. May be: temp, sal or chl"); + message(" var. Required. May be: temp, sal, chl, mld, ssh or w."); message(" time. Time moment or regular expression. If present, timeb and timee must be absent"); message(" timeb, timee. Time interval. If present, time must be absent"); message(" out. Output file. Default: out.bin"); message(" Keys for source=NEMO"); - message(" var can be only temp or sal"); + message(" var can be temp, sal, mld, ssh or w. Mld is not available for dataset=NRT6."); message(" dataset. Can be DT, NRT or NRT6. Default: DT"); - message(" layer and/or depth. Layer or depth of NEMO dataset. If depth is specified, layer is ignored. Default: layer=0"); + message(" layer and/or depth. Layer or depth of NEMO dataset. If depth is specified, layer is ignored. Both ignored, if var=mld. Default: layer=0"); message(" lonb, lone, latb, late. Required. Region of interest"); }