From 4dc5eda1dfcce74a27c2364b969d969a221d181b Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Wed, 5 Apr 2023 15:29:37 +1000 Subject: [PATCH] Merged HYCOM and NEMO classes. Configuration file name changed. URLs for data are now read from the configuration. Variables are searched for by the standard_name attribute. --- include/HYCOM.h | 392 ++----------------------------------- include/NEMO.h | 445 ++---------------------------------------- include/layereddata.h | 155 +++++++++++++++ src/ParseArgs.cpp | 2 +- src/layereddata.cpp | 290 +++++++++++++++++++++++++++ 5 files changed, 471 insertions(+), 813 deletions(-) create mode 100644 include/layereddata.h create mode 100644 src/layereddata.cpp diff --git a/include/HYCOM.h b/include/HYCOM.h index 6c2ad84..122cc97 100644 --- a/include/HYCOM.h +++ b/include/HYCOM.h @@ -1,26 +1,7 @@ #pragma once -#include "DataAdapters/ncfilealt.h" -#include "ParseArgs.h" -#include "mdatetime.h" -#include "simple2ddata.h" -#include "specfunc.h" -#include -#include -#include -#include -#include +#include "layereddata.h" -using michlib::Ceil; -using michlib::DetGeoDomain; -using michlib::errmessage; -using michlib::Floor; -using michlib::GPL; -using michlib::int2; -using michlib::MDateTime; -using michlib::NCFileA; -using michlib::ToGeoDomain; - -class HYCOMData +class HYCOMData: public LayeredData { enum Type { @@ -30,132 +11,17 @@ class HYCOMData TYPE_FORECAST }; - struct Parameters: public BaseParameters - { - size_t xb, yb, xe, ye, layer; - MString varname; - virtual ~Parameters() override = default; - }; - - std::vector nc; - std::vector urls; - std::vector depths; - std::vector times; - MString lonname, latname; - size_t nx, ny; - real lonb, latb, lone, late; - real lonstep, latstep; - Type type = TYPE_UNKNOWN; - - public: - using Data = Simple2DData; - - private: - template Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const Parameters* p) const - { - real unitmul = 1.0; - DataType fill; - real offset = 0.0, scale = 1.0; - - { - auto a_fill = f.A(name, "_FillValue"); - auto a_offset_d = f.A(name, "add_offset"); - auto a_scale_d = f.A(name, "scale_factor"); - auto a_offset_f = f.A(name, "add_offset"); - auto a_scale_f = f.A(name, "scale_factor"); - if(!a_fill) return Data(); - fill = a_fill; - if(a_offset_d) offset = a_offset_d; - if(a_scale_d) scale = a_scale_d; - if(a_offset_f) offset = a_offset_f; - if(a_scale_f) scale = a_scale_f; - } - - auto unit = f.A(name, "units"); - if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0; - - Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep); - - if(p->xb < p->xe) - { - auto var = nodepth ? f.V(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); - if(!var) return Data(); - if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data(); - - for(size_t ix = 0; ix < var.DimLen(0); ix++) - for(size_t iy = 0; iy < var.DimLen(1); iy++) - { - DataType v = var(ix, iy); - data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); - } - } - else - { - auto var1 = nodepth ? f.V(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); - auto var2 = nodepth ? f.V(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->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++) - for(size_t iy = 0; iy < var1.DimLen(1); iy++) - { - DataType v = var1(ix, iy); - data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); - } - for(size_t ix = 0; ix < var2.DimLen(0); ix++) - for(size_t iy = 0; iy < var2.DimLen(1); iy++) - { - DataType v = var2(ix, iy); - data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); - } - } - return data; - } - - static MString VName2IName(const MString& vname) - { - if(vname == "temp") return "water_temp"; - if(vname == "sal") return "salinity"; - if(vname == "ssh") return "surf_el"; - return ""; - } - - public: - Data Read(const BaseParameters* ip, size_t i) const - { - if(!isOk()) return Data(); - bool nodepth = false; - - auto p = dynamic_cast(ip); - MString name = VName2IName(p->varname); - 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, p); - if(v.Type().Id() == NC_FLOAT) return ReadVarRaw(f, name, i, nodepth, p); - } - } - return Data(); - } + Type type = TYPE_UNKNOWN; - bool HasVar(const MString vname) const + MString DataTitle() const { - MString name = VName2IName(vname); - if(name == "") return false; - - for(const auto& f: nc) + switch(type) { - auto head = f.Header(); - for(const auto& v: head.Variables()) - if(v.Name() == name) return true; + case(TYPE_REANALYSIS): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Reanalysis"; + case(TYPE_HINDCAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Hindcast"; + case(TYPE_FORECAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Forecast"; + default: return "No title"; } - return false; } public: @@ -166,251 +32,17 @@ class HYCOMData { MString dataset = args.contains("dataset") ? args.at("dataset") : "Forecast"; - nc.clear(); - + GPL.UsePrefix("HYCOM"); if(dataset == "Forecast") - { - urls = {"https://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0/FMRC/GLBy0.08_930_FMRC_best.ncd"}; type = TYPE_FORECAST; - } else if(dataset == "Hindcast") - { - urls = {"https://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0"}; type = TYPE_HINDCAST; - } else if(dataset == "Reanalysis") - { - urls.clear(); - for(size_t i = 1994; i <= 2015; i++) urls.push_back("https://tds.hycom.org/thredds/dodsC/GLBv0.08/expt_53.X/data/" + MString(i)); type = TYPE_REANALYSIS; - } else return "Unknown dataset: " + dataset; - nc.resize(urls.size()); - for(size_t i = 0; i < urls.size(); i++) - { - nc[i].Reset(urls[i] + "#cache&noprefetch"); - if(!nc[i]) - { - nc.clear(); - return "Can't connect to database"; - } - } - - auto head = nc[0].Header(); - lonname = latname = ""; - for(const auto& dim: head.Dimensions()) - { - if(dim.Name() == "lon" || dim.Name() == "longitude") - { - lonname = dim.Name(); - nx = dim.Len(); - } - if(dim.Name() == "lat" || dim.Name() == "latitude") - { - latname = dim.Name(); - ny = dim.Len(); - } - } - if(!(lonname.Exist() && latname.Exist())) - { - nc.clear(); - return "Can't find longitude/latitude"; - } - - auto rdepths = nc[0].VR("depth"); - if(!rdepths) - { - nc.clear(); - return "Can't read depths"; - } - depths.resize(rdepths.DimLen(0)); - for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i); - - auto time = nc[0].VR("time"); - if(!time) - { - nc.clear(); - return "Can't read times"; - } - - MDateTime refdate; - { - auto units = nc[0].Attribute("time", "units"); - if(!units) - { - nc.clear(); - return "Can't read refdate"; - } - MString rstr; - auto words = michlib::Split_on_words(units); - auto ci = words.begin(); - if(ci != words.end()) ci++; // skip "hours" - if(ci != words.end()) ci++; // skip "since" - if(ci != words.end()) rstr = *ci; // Day - if(ci != words.end()) ci++; - if(ci != words.end()) rstr += " " + *ci; // Hours - if(!refdate.FromString(rstr)) - { - nc.clear(); - return "Can't parse " + rstr + " to refdate"; - } - } - - times.resize(time.DimLen(0)); - for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast(time(i)) * 3600; - - auto lons = nc[0].VR(lonname); - auto lats = nc[0].VR(latname); - if(!(lons && lats)) - { - nc.clear(); - return "Can't get longitudes/latitudes"; - } - lonb = lons(0); - latb = lats(0); - lone = lons(nx - 1); - late = lats(ny - 1); - lonstep = (lone - lonb) / (nx - 1); - latstep = (late - latb) / (ny - 1); - - return ""; + SetTitle(DataTitle()); + return LayeredData::Open(args, dataset); } - - MString Info() const - { - if(!isOk()) return ""; - MString d; - for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")"; - - std::set vars; - for(const auto& f: nc) - { - auto head = f.Header(); - for(const auto& v: head.Variables()) - { - if(v.Name() == "water_temp") vars.emplace("temp"); - if(v.Name() == "salinity") vars.emplace("sal"); - if(v.Name() == "surf_el") vars.emplace("ssh"); - } - } - MString svars; - { - bool first = true; - for(const auto& v: vars) - { - svars += (first ? "" : ", ") + v; - first = false; - } - } - - // clang-format off - return - "Dataset: " + Title() + "\n" + - " Begin date: " + Time(0).ToString() + "\n" + - " End date: " + Time(NTimes()-1).ToString() + "\n" + - " Time step: " + Timestep() + " seconds\n" + - " Time moments: " + NTimes() + "\n" + - " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" + - " Grid: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" + - " Depths:" + d + "\n" + - " Supported variables: " + svars; - // clang-format on - } - - std::pair Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const - { - std::unique_ptr ppar{new struct Parameters}; - - if(!args.contains("var")) return {nullptr, "Variable not specified"}; - ppar->varname = args.at("var"); - if(!HasVar(ppar->varname)) return {nullptr, "Variable " + ppar->varname + " not exists in this dataset"}; - if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger(); - if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"}; - real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer); - - if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"}; - - { - auto dom = DetGeoDomain(lonb, lone); - real lon1 = ToGeoDomain(args.at("lonb").ToReal(), dom); - real lon2 = ToGeoDomain(args.at("lone").ToReal(), dom); - real lat1 = args.at("latb").ToReal(); - real lat2 = args.at("late").ToReal(); - - ppar->yb = static_cast(Floor((lat1 - latb) / latstep)); - ppar->ye = static_cast(Ceil((lat2 - latb) / latstep)); - if(ppar->ye > ny - 1) ppar->ye = ny - 1; - if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"}; - - ppar->xb = static_cast(Floor((lon1 - lonb) / lonstep)); - ppar->xe = static_cast(Ceil((lon2 - lonb) / lonstep)); - - if(ppar->xb == ppar->xe) return {nullptr, "Lonb must be not equal late"}; - - if(depth < 0.0 || depth > depths.back()) - ppar->layer = (depth < 0.0) ? 0 : (depths.size() - 1); - else - for(size_t i = 0; i < depths.size() - 1; i++) - { - if(depth >= depths[i] && depth <= depths[i + 1]) - { - ppar->layer = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1); - break; - } - } - } - - pars.SetParameter("variable", ppar->varname); - pars.SetParameter("depth", Depth(ppar->layer)); - pars.SetParameter("layer", ppar->layer); - pars.SetParameter("dataset", Title()); - pars.SetParameter("lonb", Lon(ppar->xb)); - pars.SetParameter("latb", Lat(ppar->yb)); - pars.SetParameter("lone", Lon(ppar->xe)); - pars.SetParameter("late", Lat(ppar->ye)); - - return {ppar.release(), ""}; - } - - MString Dump(const struct Parameters* ppar) const - { - // clang-format off - return - "Current settings:\n" + MString() + - " Longitudes: from " + Lon(ppar->xb) + " (" + ppar->xb + ") to "+ Lon(ppar->xe) + " (" + ppar->xe + ")\n" + - " Latitudes: from " + Lat(ppar->yb) + " (" + ppar->yb + ") to "+ Lat(ppar->ye) + " (" + ppar->ye + ")\n" + - " Depth: layer " + ppar->layer + ", depth " + Depth(ppar->layer) + " m\n"; - // clang-format on - } - - 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; } - real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; } - time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; } - MDateTime Time(size_t i) const - { - if(!isOk() || i >= times.size()) return MDateTime(); - return times[i]; - } - - 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(); } - - MString Title() const - { - switch(type) - { - case(TYPE_REANALYSIS): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Reanalysis"; - case(TYPE_HINDCAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Hindcast"; - case(TYPE_FORECAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Forecast"; - default: return "No title"; - } - } - - static real Fillval() { return Data::Fillval(); } }; diff --git a/include/NEMO.h b/include/NEMO.h index f63578a..d2b4d55 100644 --- a/include/NEMO.h +++ b/include/NEMO.h @@ -1,26 +1,7 @@ #pragma once -#include "DataAdapters/ncfilealt.h" -#include "ParseArgs.h" -#include "mdatetime.h" -#include "simple2ddata.h" -#include "specfunc.h" -#include -#include -#include -#include -#include +#include "layereddata.h" -using michlib::Ceil; -using michlib::DetGeoDomain; -using michlib::errmessage; -using michlib::Floor; -using michlib::GPL; -using michlib::int2; -using michlib::MDateTime; -using michlib::NCFileA; -using michlib::ToGeoDomain; - -class NEMOData +class NEMOData: public LayeredData { enum Type { @@ -30,170 +11,17 @@ class NEMOData TYPE_NRT6 }; - struct Parameters: public BaseParameters - { - size_t xb, yb, xe, ye, layer; - MString varname; - virtual ~Parameters() override = default; - }; - - std::vector nc; - std::vector urls; - std::vector depths; - std::vector times; - MString lonname, latname; - size_t nx, ny; - real lonb, latb, lone, late; - real lonstep, latstep; - Type type = TYPE_UNKNOWN; - - class EnvVar - { - MString name, oldvalue; - bool activated, saved; - - public: - EnvVar(): activated(false) {} - ~EnvVar() { Deactivate(); } - - void Activate(const MString& var, const MString& val) - { - if(activated) Deactivate(); - name = var; - char* curval = getenv(name.Buf()); - if(nullptr == curval) - saved = false; - else - { - oldvalue = curval; - saved = true; - } - setenv(name.Buf(), val.Buf(), 1); - } - void Deactivate() - { - if(!activated) return; - if(saved) - setenv(name.Buf(), oldvalue.Buf(), 1); - else - unsetenv(name.Buf()); - activated = false; - } - }; - - EnvVar proxy; - - public: - using Data = Simple2DData; + Type type = TYPE_UNKNOWN; - private: - template Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const Parameters* p) const + MString DataTitle() const { - real unitmul = 1.0; - DataType fill; - real offset = 0.0, scale = 1.0; - - { - auto a_fill = f.A(name, "_FillValue"); - auto a_offset_d = f.A(name, "add_offset"); - auto a_scale_d = f.A(name, "scale_factor"); - auto a_offset_f = f.A(name, "add_offset"); - auto a_scale_f = f.A(name, "scale_factor"); - if(!a_fill) return Data(); - fill = a_fill; - if(a_offset_d) offset = a_offset_d; - if(a_scale_d) scale = a_scale_d; - if(a_offset_f) offset = a_offset_f; - if(a_scale_f) scale = a_scale_f; - } - - auto unit = f.A(name, "units"); - if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0; - - Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep); - - if(p->xb < p->xe) - { - auto var = nodepth ? f.V(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); - if(!var) return Data(); - if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data(); - - for(size_t ix = 0; ix < var.DimLen(0); ix++) - for(size_t iy = 0; iy < var.DimLen(1); iy++) - { - DataType v = var(ix, iy); - data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); - } - } - else - { - auto var1 = nodepth ? f.V(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); - auto var2 = nodepth ? f.V(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->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++) - for(size_t iy = 0; iy < var1.DimLen(1); iy++) - { - DataType v = var1(ix, iy); - data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); - } - for(size_t ix = 0; ix < var2.DimLen(0); ix++) - for(size_t iy = 0; iy < var2.DimLen(1); iy++) - { - DataType v = var2(ix, iy); - data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); - } - } - return data; - } - - static MString VName2IName(const MString& vname) - { - if(vname == "temp") return "thetao"; - if(vname == "sal") return "so"; - if(vname == "mld") return "mlotst"; - if(vname == "ssh") return "zos"; - if(vname == "w") return "wo"; - return ""; - } - - public: - Data Read(const BaseParameters* ip, size_t i) const - { - if(!isOk()) return Data(); - bool nodepth = false; - - auto p = dynamic_cast(ip); - MString name = VName2IName(p->varname); - 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, p); - if(v.Type().Id() == NC_FLOAT) return ReadVarRaw(f, name, i, nodepth, p); - } - } - return Data(); - } - - bool HasVar(const MString vname) const - { - MString name = VName2IName(vname); - if(name == "") return false; - - for(const auto& f: nc) + switch(type) { - auto head = f.Header(); - for(const auto& v: head.Variables()) - if(v.Name() == name) return true; + case(TYPE_DT): return "NEMO Delayed time, daily mean (DT)"; + case(TYPE_NRT): return "NEMO Near-real time, daily mean (NRT)"; + case(TYPE_NRT6): return "NEMO Near-real time, 6h resolution (NRT6)"; + default: return "No title"; } - return false; } public: @@ -202,266 +30,19 @@ class NEMOData // TODO: RetVal MString Open(const CLArgs& args) { - MString oldprefix = GPL.UsePrefix("AVISO"); - MString dataset = args.contains("dataset") ? args.at("dataset") : "DT"; - MString cred = GPL.ParameterSValue("COPERNICUS_USER", ""); - MString proxyurl = GPL.ParameterSValue("COPERNICUS_PROXY", ""); - - GPL.UsePrefix(oldprefix); - - nc.clear(); - if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl); + MString dataset = args.contains("dataset") ? args.at("dataset") : "DT"; + GPL.UsePrefix("NEMO"); if(dataset == "DT") - { - urls = {"https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m"}; type = TYPE_DT; - } else if(dataset == "NRT") - { - 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"}; type = TYPE_NRT; - } else if(dataset == "NRT6") - { - 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"}; type = TYPE_NRT6; - } else return "Unknown dataset: " + dataset; - nc.resize(urls.size()); - for(size_t i = 0; i < urls.size(); i++) - { - nc[i].Reset(urls[i] + "#cache&noprefetch"); - if(!nc[i]) - { - nc.clear(); - return "Can't connect to database"; - } - } - - auto head = nc[0].Header(); - lonname = latname = ""; - for(const auto& dim: head.Dimensions()) - { - if(dim.Name() == "lon" || dim.Name() == "longitude") - { - lonname = dim.Name(); - nx = dim.Len(); - } - if(dim.Name() == "lat" || dim.Name() == "latitude") - { - latname = dim.Name(); - ny = dim.Len(); - } - } - if(!(lonname.Exist() && latname.Exist())) - { - nc.clear(); - return "Can't find longitude/latitude"; - } - - auto rdepths = nc[0].VR("depth"); - if(!rdepths) - { - nc.clear(); - return "Can't read depths"; - } - depths.resize(rdepths.DimLen(0)); - for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i); - - auto time = nc[0].VR("time"); - if(!time) - { - nc.clear(); - return "Can't read times"; - } - - MDateTime refdate; - { - auto units = nc[0].Attribute("time", "units"); - if(!units) - { - nc.clear(); - return "Can't read refdate"; - } - MString rstr; - auto words = michlib::Split_on_words(units); - auto ci = words.begin(); - if(ci != words.end()) ci++; // skip "hours" - if(ci != words.end()) ci++; // skip "since" - if(ci != words.end()) rstr = *ci; // Day - if(ci != words.end()) ci++; - if(ci != words.end()) rstr += " " + *ci; // Hours - if(!refdate.FromString(rstr)) - { - nc.clear(); - return "Can't parse " + rstr + " to refdate"; - } - } - - times.resize(time.DimLen(0)); - for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast(time(i)) * 3600; - - auto lons = nc[0].VR(lonname); - auto lats = nc[0].VR(latname); - if(!(lons && lats)) - { - nc.clear(); - return "Can't get longitudes/latitudes"; - } - lonb = lons(0); - latb = lats(0); - lone = lons(nx - 1); - late = lats(ny - 1); - lonstep = (lone - lonb) / (nx - 1); - latstep = (late - latb) / (ny - 1); - - return ""; + SetTitle(DataTitle()); + return LayeredData::Open(args, dataset); } - - MString Info() const - { - if(!isOk()) return ""; - MString d; - for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")"; - - std::set vars; - for(const auto& f: nc) - { - auto head = f.Header(); - for(const auto& v: head.Variables()) - { - if(v.Name() == "thetao") vars.emplace("temp"); - if(v.Name() == "so") vars.emplace("sal"); - if(v.Name() == "mlotst") vars.emplace("mld"); - if(v.Name() == "zos") vars.emplace("ssh"); - if(v.Name() == "wo") vars.emplace("w"); - } - } - MString svars; - { - bool first = true; - for(const auto& v: vars) - { - svars += (first ? "" : ", ") + v; - first = false; - } - } - - // clang-format off - return - "Dataset: " + Title() + "\n" + - " Begin date: " + Time(0).ToString() + "\n" + - " End date: " + Time(NTimes()-1).ToString() + "\n" + - " Time step: " + Timestep() + " seconds\n" + - " Time moments: " + NTimes() + "\n" + - " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" + - " Grid: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" + - " Depths:" + d + "\n" + - " Supported variables: " + svars; - // clang-format on - } - - std::pair Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const - { - std::unique_ptr ppar{new struct Parameters}; - - if(!args.contains("var")) return {nullptr, "Variable not specified"}; - ppar->varname = args.at("var"); - if(!HasVar(ppar->varname)) return {nullptr, "Variable " + ppar->varname + " not exists in this dataset"}; - if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger(); - if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"}; - real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer); - - if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"}; - - { - auto dom = DetGeoDomain(lonb, lone); - real lon1 = ToGeoDomain(args.at("lonb").ToReal(), dom); - real lon2 = ToGeoDomain(args.at("lone").ToReal(), dom); - real lat1 = args.at("latb").ToReal(); - real lat2 = args.at("late").ToReal(); - - ppar->yb = static_cast(Floor((lat1 - latb) / latstep)); - ppar->ye = static_cast(Ceil((lat2 - latb) / latstep)); - if(ppar->ye > ny - 1) ppar->ye = ny - 1; - if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"}; - - ppar->xb = static_cast(Floor((lon1 - lonb) / lonstep)); - ppar->xe = static_cast(Ceil((lon2 - lonb) / lonstep)); - - if(ppar->xb == ppar->xe) return {nullptr, "Lonb must be not equal late"}; - - if(depth < 0.0 || depth > depths.back()) - ppar->layer = (depth < 0.0) ? 0 : (depths.size() - 1); - else - for(size_t i = 0; i < depths.size() - 1; i++) - { - if(depth >= depths[i] && depth <= depths[i + 1]) - { - ppar->layer = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1); - break; - } - } - } - - pars.SetParameter("variable", ppar->varname); - pars.SetParameter("depth", Depth(ppar->layer)); - pars.SetParameter("layer", ppar->layer); - pars.SetParameter("dataset", Title()); - pars.SetParameter("lonb", Lon(ppar->xb)); - pars.SetParameter("latb", Lat(ppar->yb)); - pars.SetParameter("lone", Lon(ppar->xe)); - pars.SetParameter("late", Lat(ppar->ye)); - - return {ppar.release(), ""}; - } - - MString Dump(const struct Parameters* ppar) const - { - // clang-format off - return - "Current settings:\n" + MString() + - " Longitudes: from " + Lon(ppar->xb) + " (" + ppar->xb + ") to "+ Lon(ppar->xe) + " (" + ppar->xe + ")\n" + - " Latitudes: from " + Lat(ppar->yb) + " (" + ppar->yb + ") to "+ Lat(ppar->ye) + " (" + ppar->ye + ")\n" + - " Depth: layer " + ppar->layer + ", depth " + Depth(ppar->layer) + " m\n"; - // clang-format on - } - - 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; } - real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; } - time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; } - MDateTime Time(size_t i) const - { - if(!isOk() || i >= times.size()) return MDateTime(); - return times[i]; - } - - 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(); } - - MString Title() const - { - switch(type) - { - case(TYPE_DT): return "NEMO Delayed time, daily mean (DT)"; - case(TYPE_NRT): return "NEMO Near-real time, daily mean (NRT)"; - case(TYPE_NRT6): return "NEMO Near-real time, 6h resolution (NRT6)"; - default: return "No title"; - } - } - - static real Fillval() { return Data::Fillval(); } }; diff --git a/include/layereddata.h b/include/layereddata.h new file mode 100644 index 0000000..ab72d2f --- /dev/null +++ b/include/layereddata.h @@ -0,0 +1,155 @@ +#pragma once +#include "DataAdapters/ncfilealt.h" +#include "ParseArgs.h" +#include "mdatetime.h" +#include "simple2ddata.h" +#include +#include + +using michlib::Ceil; +using michlib::DetGeoDomain; +using michlib::Floor; +using michlib::GPL; +using michlib::int2; +using michlib::MDateTime; +using michlib::MString; +using michlib::NCFileA; +using michlib::ToGeoDomain; + +class LayeredData +{ + public: + using Data = Simple2DData; + + private: + std::vector nc; + std::vector urls; + std::vector depths; + std::vector times; + MString lonname, latname; + size_t nx, ny; + real lonb, latb, lone, late; + real lonstep, latstep; + MString title; + + class EnvVar + { + MString name, oldvalue; + bool activated, saved; + + public: + EnvVar(): activated(false) {} + ~EnvVar() { Deactivate(); } + + void Activate(const MString& var, const MString& val) + { + if(activated) Deactivate(); + name = var; + char* curval = getenv(name.Buf()); + if(nullptr == curval) + saved = false; + else + { + oldvalue = curval; + saved = true; + } + setenv(name.Buf(), val.Buf(), 1); + } + void Deactivate() + { + if(!activated) return; + if(saved) + setenv(name.Buf(), oldvalue.Buf(), 1); + else + unsetenv(name.Buf()); + activated = false; + } + }; + + EnvVar proxy; + + protected: + struct Parameters: public BaseParameters + { + size_t xb, yb, xe, ye, layer; + MString varname; + virtual ~Parameters() override = default; + }; + + // TODO: RetVal + MString Open(const CLArgs& args, const MString& dataset); + + void SetTitle(const MString& newtitle) { title = newtitle; } + + public: + MString Info() const; + + std::pair Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const; + + Data Read(const BaseParameters* ip, size_t i) const; + + bool isOk() const { return nc.size() > 0; } + + explicit operator bool() const { return nc.size() > 0; } + + real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; } + + 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; } + + size_t NDepths() const { return depths.size(); } + + size_t NTimes() const { return times.size(); } + + MDateTime Time(size_t i) const + { + if(!isOk() || i >= times.size()) return MDateTime(); + return times[i]; + } + + time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; } + + MString Title() const { return title; } + + MString Dump(const struct Parameters* ppar) const + { + // clang-format off + return + "Current settings:\n" + MString() + + " Longitudes: from " + Lon(ppar->xb) + " (" + ppar->xb + ") to "+ Lon(ppar->xe) + " (" + ppar->xe + ")\n" + + " Latitudes: from " + Lat(ppar->yb) + " (" + ppar->yb + ") to "+ Lat(ppar->ye) + " (" + ppar->ye + ")\n" + + " Depth: layer " + ppar->layer + ", depth " + Depth(ppar->layer) + " m\n"; + // clang-format on + } + + private: + template Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const struct Parameters* p) const; + + std::pair VarNameLoc(const MString vname) const + { + for(size_t i = 0; i < nc.size(); i++) + { + auto head = nc[i].Header(); + for(const auto& v: head.Variables()) + { + auto stname = nc[i].A(v.Name(), "standard_name"); + if(!stname) continue; + if(StName2Name(stname) == vname) return {v.Name(), i}; + } + } + return {"", 0}; + } + + static MString StName2Name(const MString& stname) + { + if(stname == "sea_water_potential_temperature") return "temp"; + if(stname == "sea_water_temperature") return "temp"; + if(stname == "sea_water_salinity") return "sal"; + if(stname == "ocean_mixed_layer_thickness_defined_by_sigma_theta") return "mld"; + if(stname == "sea_surface_height_above_geoid") return "ssh"; + if(stname == "sea_surface_elevation") return "ssh"; + if(stname == "upward_sea_water_velocity") return "w"; + return ""; + } +}; diff --git a/src/ParseArgs.cpp b/src/ParseArgs.cpp index 517fc86..9b2d623 100644 --- a/src/ParseArgs.cpp +++ b/src/ParseArgs.cpp @@ -28,7 +28,7 @@ CLArgs ParseArgs(int argc, char** argv) } SList sl; - michlib_internal::ParseParameterFile("/etc/michlibconf", sl, 0); + michlib_internal::ParseParameterFile("/etc/odm.conf", sl, 0); return out; } diff --git a/src/layereddata.cpp b/src/layereddata.cpp new file mode 100644 index 0000000..35d31d3 --- /dev/null +++ b/src/layereddata.cpp @@ -0,0 +1,290 @@ +#define MICHLIB_NOSOURCE +#include "layereddata.h" + +MString LayeredData::Info() const +{ + if(!isOk()) return ""; + MString d; + for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")"; + + std::set vars; + for(const auto& f: nc) + { + auto head = f.Header(); + for(const auto& v: head.Variables()) + { + auto ret = f.A(v.Name(), "standard_name"); + if(!ret) continue; + vars.emplace(StName2Name(ret)); + } + } + MString svars; + { + bool first = true; + for(const auto& v: vars) + { + svars += (first ? "" : ", ") + v; + first = false; + } + } + + // clang-format off + return + "Dataset: " + Title() + "\n" + + " Begin date: " + Time(0).ToString() + "\n" + + " End date: " + Time(NTimes()-1).ToString() + "\n" + + " Time step: " + Timestep() + " seconds\n" + + " Time moments: " + NTimes() + "\n" + + " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" + + " Grid: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" + + " Depths:" + d + "\n" + + " Supported variables: " + svars; + // clang-format on +} + +MString LayeredData::Open(const CLArgs& args, const MString& dataset) +{ + nc.clear(); + MString proxyurl = GPL.ParameterSValue("USEPROXY", ""); + if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl); + + urls.clear(); + size_t i = 1; + while(true) + { + MString url = GPL.ParameterSValue(dataset + "_URL" + i, ""); + if(url.Exist()) + urls.emplace_back(std::move(url)); + else + break; + i++; + } + if(urls.size() == 0) return "No urls for dataset " + dataset + " specified in config"; + + nc.resize(urls.size()); + for(size_t i = 0; i < urls.size(); i++) + { + nc[i].Reset(urls[i] + "#cache&noprefetch"); + if(!nc[i]) + { + nc.clear(); + return "Can't connect to database"; + } + } + + auto head = nc[0].Header(); + lonname = latname = ""; + for(const auto& dim: head.Dimensions()) + { + if(dim.Name() == "lon" || dim.Name() == "longitude") + { + lonname = dim.Name(); + nx = dim.Len(); + } + if(dim.Name() == "lat" || dim.Name() == "latitude") + { + latname = dim.Name(); + ny = dim.Len(); + } + } + if(!(lonname.Exist() && latname.Exist())) + { + nc.clear(); + return "Can't find longitude/latitude"; + } + + auto rdepths = nc[0].VR("depth"); + if(!rdepths) + { + nc.clear(); + return "Can't read depths"; + } + depths.resize(rdepths.DimLen(0)); + for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i); + + auto time = nc[0].VR("time"); + if(!time) + { + nc.clear(); + return "Can't read times"; + } + + MDateTime refdate; + { + auto units = nc[0].Attribute("time", "units"); + if(!units) + { + nc.clear(); + return "Can't read refdate"; + } + MString rstr; + auto words = michlib::Split_on_words(units); + auto ci = words.begin(); + if(ci != words.end()) ci++; // skip "hours" + if(ci != words.end()) ci++; // skip "since" + if(ci != words.end()) rstr = *ci; // Day + if(ci != words.end()) ci++; + if(ci != words.end()) rstr += " " + *ci; // Hours + if(!refdate.FromString(rstr)) + { + nc.clear(); + return "Can't parse " + rstr + " to refdate"; + } + } + + times.resize(time.DimLen(0)); + for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast(time(i)) * 3600; + + auto lons = nc[0].VR(lonname); + auto lats = nc[0].VR(latname); + if(!(lons && lats)) + { + nc.clear(); + return "Can't get longitudes/latitudes"; + } + lonb = lons(0); + latb = lats(0); + lone = lons(nx - 1); + late = lats(ny - 1); + lonstep = (lone - lonb) / (nx - 1); + latstep = (late - latb) / (ny - 1); + + return ""; +} + +std::pair LayeredData::Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const +{ + std::unique_ptr ppar{new struct Parameters}; + + if(!args.contains("var")) return {nullptr, "Variable not specified"}; + ppar->varname = args.at("var"); + if(!VarNameLoc(ppar->varname).first.Exist()) return {nullptr, "Variable " + ppar->varname + " not exists in this dataset"}; + if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger(); + if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"}; + real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer); + + if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"}; + + { + auto dom = DetGeoDomain(lonb, lone); + real lon1 = ToGeoDomain(args.at("lonb").ToReal(), dom); + real lon2 = ToGeoDomain(args.at("lone").ToReal(), dom); + real lat1 = args.at("latb").ToReal(); + real lat2 = args.at("late").ToReal(); + + ppar->yb = static_cast(Floor((lat1 - latb) / latstep)); + ppar->ye = static_cast(Ceil((lat2 - latb) / latstep)); + if(ppar->ye > ny - 1) ppar->ye = ny - 1; + if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"}; + + ppar->xb = static_cast(Floor((lon1 - lonb) / lonstep)); + ppar->xe = static_cast(Ceil((lon2 - lonb) / lonstep)); + + if(ppar->xb == ppar->xe) return {nullptr, "Lonb must be not equal late"}; + + if(depth < 0.0 || depth > depths.back()) + ppar->layer = (depth < 0.0) ? 0 : (depths.size() - 1); + else + for(size_t i = 0; i < depths.size() - 1; i++) + { + if(depth >= depths[i] && depth <= depths[i + 1]) + { + ppar->layer = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1); + break; + } + } + } + + pars.SetParameter("variable", ppar->varname); + pars.SetParameter("depth", Depth(ppar->layer)); + pars.SetParameter("layer", ppar->layer); + pars.SetParameter("dataset", Title()); + pars.SetParameter("lonb", Lon(ppar->xb)); + pars.SetParameter("latb", Lat(ppar->yb)); + pars.SetParameter("lone", Lon(ppar->xe)); + pars.SetParameter("late", Lat(ppar->ye)); + + return {ppar.release(), ""}; +} + +LayeredData::Data LayeredData::Read(const BaseParameters* ip, size_t i) const +{ + if(!isOk()) return Data(); + bool nodepth = false; + + auto p = dynamic_cast(ip); + auto [name, id] = VarNameLoc(p->varname); + auto head = nc[id].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(nc[id], name, i, nodepth, p); + if(v.Type().Id() == NC_FLOAT) return ReadVarRaw(nc[id], name, i, nodepth, p); + } + + return Data(); +} + +template LayeredData::Data LayeredData::ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const struct LayeredData::Parameters* p) const +{ + real unitmul = 1.0; + DataType fill; + real offset = 0.0, scale = 1.0; + + { + auto a_fill = f.A(name, "_FillValue"); + auto a_offset_d = f.A(name, "add_offset"); + auto a_scale_d = f.A(name, "scale_factor"); + auto a_offset_f = f.A(name, "add_offset"); + auto a_scale_f = f.A(name, "scale_factor"); + if(!a_fill) return Data(); + fill = a_fill; + if(a_offset_d) offset = a_offset_d; + if(a_scale_d) scale = a_scale_d; + if(a_offset_f) offset = a_offset_f; + if(a_scale_f) scale = a_scale_f; + } + + auto unit = f.A(name, "units"); + if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0; + + Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep); + + if(p->xb < p->xe) + { + auto var = nodepth ? f.V(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); + if(!var) return Data(); + if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data(); + + for(size_t ix = 0; ix < var.DimLen(0); ix++) + for(size_t iy = 0; iy < var.DimLen(1); iy++) + { + DataType v = var(ix, iy); + data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); + } + } + else + { + auto var1 = nodepth ? f.V(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); + auto var2 = nodepth ? f.V(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->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++) + for(size_t iy = 0; iy < var1.DimLen(1); iy++) + { + DataType v = var1(ix, iy); + data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); + } + for(size_t ix = 0; ix < var2.DimLen(0); ix++) + for(size_t iy = 0; iy < var2.DimLen(1); iy++) + { + DataType v = var2(ix, iy); + data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); + } + } + return data; +}