From 5359665083deef7e7ed130459ff0c47f0b58066f Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Fri, 7 Apr 2023 15:48:10 +1000 Subject: [PATCH] Added support for time-splitted data (HYCOM Reanalysis) --- include/layereddata.h | 93 ++++++++++++++++++++++++++++--- src/layereddata.cpp | 124 +++++++++++++++++------------------------- 2 files changed, 136 insertions(+), 81 deletions(-) diff --git a/include/layereddata.h b/include/layereddata.h index a69d97e..d75eb91 100644 --- a/include/layereddata.h +++ b/include/layereddata.h @@ -3,6 +3,7 @@ #include "ParseArgs.h" #include "mdatetime.h" #include "simple2ddata.h" +#include #include #include @@ -22,8 +23,67 @@ class LayeredData using Data = Simple2DData; private: - std::vector nc; - std::vector urls; + class NC + { + MString url; + NCFileA nc; + std::vector times; + + public: + NC(MString&& newurl): url(std::move(newurl)) { nc.Reset(url + "#cache&noprefetch"); } + + MString ReadTimes() + { + if(!nc) return "File not open"; + auto time = nc.VR("time"); + if(!time) return "Can't read times"; + + MDateTime refdate; + { + auto units = nc.Attribute("time", "units"); + if(!units) 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)) 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; + return ""; + } + + const NCFileA* operator->() const { return &nc; } + explicit operator bool() const { return nc; } + MDateTime Begin() const { return times.front(); } + MDateTime End() const { return times.back(); } + + const std::vector& Times() const { return times; } + size_t Index(MDateTime tm) const + { + if(tm < Begin() || tm > End()) return 0; + size_t b = 0, e = times.size() - 1; + if(tm == times[b]) return b + 1; + if(tm == times[e]) return e + 1; + while(e - b > 1) + { + size_t c = (e + b) / 2; + if(tm == times[c]) return c + 1; + if(tm > times[c]) + b = c; + else + e = c; + } + return 0; + } + }; + std::vector nc; std::vector depths; std::vector times; MString lonname, latname; @@ -124,21 +184,38 @@ class LayeredData } private: - template Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const struct Parameters* p) const; + template Data ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct Parameters* p) const; + + bool HaveVar(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 true; + } + } + return false; + } - std::pair VarNameLoc(const MString vname) const + std::tuple VarNameLoc(const MString vname, MDateTime tm) const { for(size_t i = 0; i < nc.size(); i++) { - auto head = nc[i].Header(); + auto head = nc[i]->Header(); + auto tind = nc[i].Index(tm); + if(tind == 0) continue; for(const auto& v: head.Variables()) { - auto stname = nc[i].A(v.Name(), "standard_name"); + auto stname = nc[i]->A(v.Name(), "standard_name"); if(!stname) continue; - if(StName2Name(stname) == vname) return {v.Name(), i}; + if(StName2Name(stname) == vname) return {v.Name(), i, tind - 1}; } } - return {"", 0}; + return {"", 0, 0}; } static MString StName2Name(const MString& stname) diff --git a/src/layereddata.cpp b/src/layereddata.cpp index 3dea72a..3ff2097 100644 --- a/src/layereddata.cpp +++ b/src/layereddata.cpp @@ -10,12 +10,12 @@ MString LayeredData::Info() const std::set vars; for(const auto& f: nc) { - auto head = f.Header(); + auto head = f->Header(); for(const auto& v: head.Variables()) { - auto ret = f.A(v.Name(), "standard_name"); + auto ret = f->A(v.Name(), "standard_name"); if(!ret) continue; - vars.emplace(StName2Name(ret)); + if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret)); } } MString svars; @@ -48,31 +48,34 @@ MString LayeredData::Open(const MString& dataset) MString proxyurl = GPL.ParameterSValue("USEPROXY", ""); if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl); - urls.clear(); + nc.clear(); size_t i = 1; while(true) { MString url = GPL.ParameterSValue(dataset + "_URL" + i, ""); if(url.Exist()) - urls.emplace_back(std::move(url)); + { + //michlib::message("Open "+url); + nc.emplace_back(std::move(url)); + if(!nc.back()) + { + nc.clear(); + return "Can't connect to url " + url; + } + MString ret = nc.back().ReadTimes(); + if(ret.Exist()) + { + nc.clear(); + return "Can't connect to url " + url + ": " + ret; + } + } else break; i++; } - if(urls.size() == 0) return "No urls for dataset " + dataset + " specified in config"; + if(nc.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(); + auto head = nc[0]->Header(); lonname = latname = ""; for(const auto& dim: head.Dimensions()) { @@ -93,7 +96,15 @@ MString LayeredData::Open(const MString& dataset) return "Can't find longitude/latitude"; } - auto rdepths = nc[0].VR("depth"); + // Read times + { + for(const auto& f: nc) times.insert(times.end(), f.Times().begin(), f.Times().end()); + std::sort(times.begin(), times.end()); + auto last = std::unique(times.begin(), times.end()); + times.erase(last, times.end()); + } + + auto rdepths = nc[0]->VR("depth"); if(!rdepths) { nc.clear(); @@ -102,41 +113,8 @@ MString LayeredData::Open(const MString& dataset) 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); + auto lons = nc[0]->VR(lonname); + auto lats = nc[0]->VR(latname); if(!(lons && lats)) { nc.clear(); @@ -158,7 +136,7 @@ std::pair LayeredData::Parameters(michlib_intern 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(!HaveVar(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); @@ -212,32 +190,32 @@ 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(); + auto p = dynamic_cast(ip); + auto [name, id, tid] = VarNameLoc(p->varname, times[i]); + 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); + if(v.Type().Id() == NC_SHORT) return ReadVarRaw(nc[id], name, tid, nodepth, p); + if(v.Type().Id() == NC_FLOAT) return ReadVarRaw(nc[id], name, tid, 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 +template LayeredData::Data LayeredData::ReadVarRaw(const NC& 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"); + 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; @@ -246,15 +224,15 @@ template LayeredData::Data LayeredData::ReadVarRaw(const NCFileA if(a_scale_f) scale = a_scale_f; } - auto unit = f.A(name, "units"); + 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}); + 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(); @@ -267,10 +245,10 @@ template LayeredData::Data LayeredData::ReadVarRaw(const NCFileA } 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}); + 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++)