#pragma once #include "DataAdapters/ncfilealt.h" #include "ParseArgs.h" #include "mdatetime.h" #include "simple2ddata.h" #include "specfunc.h" #include #include #include #include #include using michlib::Ceil; using michlib::errmessage; using michlib::Floor; using michlib::GPL; using michlib::int2; using michlib::MDateTime; using michlib::NCFileA; using michlib::ToGeoDomainNeg; class NEMOData { enum Type { TYPE_UNKNOWN, TYPE_DT, TYPE_NRT, TYPE_NRT6 }; struct Parameters: public BaseParameters { size_t xb, yb, xe, ye, layer; MString varname; virtual ~Parameters() override = default; }; std::vector nc; std::vector depths; std::vector times; 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; 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 = 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" || 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, {"longitude", p->xb, p->xe - p->xb + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) : f.V(name, {"longitude", p->xb, p->xe - p->xb + 1}, {"latitude", 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, {"longitude", p->xb}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) : f.V(name, {"longitude", p->xb}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); auto var2 = nodepth ? f.V(name, {"longitude", 0, p->xe + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) : f.V(name, {"longitude", 0, p->xe + 1}, {"latitude", 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) { auto head = f.Header(); for(const auto& v: head.Variables()) if(v.Name() == name) return true; } return false; } public: NEMOData() = default; // 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); std::vector urls; 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]); if(!nc[i]) { nc.clear(); return "Can't connect to database"; } } auto rdepths = nc[0].V("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 timeD = nc[0].V("time"); auto timeF = nc[0].V("time"); if(!(timeD || timeF)) { nc.clear(); return "Can't read times"; } MDateTime refdate("1950-01-01"); timeD ? times.resize(timeD.DimLen(0)) : times.resize(timeF.DimLen(0)); for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast(timeD ? timeD(i) : timeF(i)) * 3600; auto nlon = nc[0].D("longitude"); auto nlat = nc[0].D("latitude"); if(!(nlon && nlat)) { nc.clear(); return "Can't get number of longitudes/latitudes"; } nx = nlon; ny = nlat; auto lons = nc[0].V("longitude"); auto lats = nc[0].V("latitude"); 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 ""; } 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)"}; { real lon1 = ToGeoDomainNeg(args.at("lonb").ToReal()); real lon2 = ToGeoDomainNeg(args.at("lone").ToReal()); 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(); } };