From 0db91e94c3248c3edb9b8e96db760c92718d8725 Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Thu, 23 Mar 2023 16:08:47 +1000 Subject: [PATCH] HYCOM support (dataset=analysis not work for now) --- include/HYCOM.h | 404 ++++++++++++++++++++++++++++++++++++++++++++++++ include/odm.h | 10 +- 2 files changed, 413 insertions(+), 1 deletion(-) create mode 100644 include/HYCOM.h diff --git a/include/HYCOM.h b/include/HYCOM.h new file mode 100644 index 0000000..9ca3984 --- /dev/null +++ b/include/HYCOM.h @@ -0,0 +1,404 @@ +#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::ToGeoDomainPos; + +class HYCOMData +{ + enum Type + { + TYPE_UNKNOWN, + TYPE_REANALYSIS, + TYPE_HINDCAST, + 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; + 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, {"lon", p->xb, p->xe - p->xb + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(name, {"lon", p->xb, p->xe - p->xb + 1}, {"lat", 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, {"lon", p->xb}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(name, {"lon", p->xb}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); + auto var2 = nodepth ? f.V(name, {"lon", 0, p->xe + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(name, {"lon", 0, p->xe + 1}, {"lat", 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(); + } + + 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: + HYCOMData() = default; + + // TODO: RetVal + MString Open(const CLArgs& args) + { + MString dataset = args.contains("dataset") ? args.at("dataset") : "Forecast"; + + nc.clear(); + + 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 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; + { + 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"; + } + } + + 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("lon"); + auto nlat = nc[0].D("lat"); + if(!(nlon && nlat)) + { + nc.clear(); + return "Can't get number of longitudes/latitudes"; + } + nx = nlon; + ny = nlat; + + auto lons = nc[0].V("lon"); + auto lats = nc[0].V("lat"); + 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() == "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)"}; + + { + real lon1 = ToGeoDomainPos(args.at("lonb").ToReal()); + real lon2 = ToGeoDomainPos(args.at("lone").ToReal()); + real lat1 = ToGeoDomainPos(args.at("latb").ToReal()); + real lat2 = ToGeoDomainPos(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/odm.h b/include/odm.h index 79e7b98..50a4e70 100644 --- a/include/odm.h +++ b/include/odm.h @@ -1,6 +1,7 @@ #pragma once #include "BFileW.h" #include "NEMO.h" +#include "HYCOM.h" #include #include @@ -9,7 +10,7 @@ using michlib::errmessage; using michlib::GPL; using michlib::message; -using DataVariants = std::variant; +using DataVariants = std::variant; template concept InfoSupported = requires(T t) { @@ -143,6 +144,13 @@ class Data: public DataVariants if(res.Exist()) return "Can't open source " + src + ":\n" + res; *this = Data(std::move(data)); } + else if(src == "HYCOM") + { + HYCOMData data; + auto res = data.Open(args); + if(res.Exist()) return "Can't open source " + src + ":\n" + res; + *this = Data(std::move(data)); + } else return "Unknown source: " + src; return "";