From 712be944663b5cff20e22280a02c9bed9e7acd8b Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Wed, 19 Apr 2023 15:38:33 +1000 Subject: [PATCH] Added source AVISOLOCAL and made a some code rearrangement. --- include/AVISOLOCAL.h | 71 ++++++++++++ include/layereddata.h | 72 +++--------- include/ncfuncs.h | 40 +++++++ include/odm.h | 10 +- include/simple2ddata.h | 3 + src/AVISOLOCAL.cpp | 258 +++++++++++++++++++++++++++++++++++++++++ src/layereddata.cpp | 125 ++++++++++---------- src/ncfuncs.cpp | 148 +++++++++++++++++++++++ src/odm.cpp | 11 +- 9 files changed, 609 insertions(+), 129 deletions(-) create mode 100644 include/AVISOLOCAL.h create mode 100644 include/ncfuncs.h create mode 100644 src/AVISOLOCAL.cpp create mode 100644 src/ncfuncs.cpp diff --git a/include/AVISOLOCAL.h b/include/AVISOLOCAL.h new file mode 100644 index 0000000..9948aab --- /dev/null +++ b/include/AVISOLOCAL.h @@ -0,0 +1,71 @@ +#pragma once +#include "GPL.h" +#include "ParseArgs.h" +#include "mregex.h" +#include "ncfuncs.h" +#include "simple2ddata.h" +#include "uvdata.h" +#include +#include + +using michlib::Ceil; +using michlib::DetGeoDomain; +using michlib::Floor; +using michlib::GPL; +using michlib::int2; +using michlib::int4; +using michlib::MString; +using michlib::real; +using michlib::RegExp; + +class AVISOLOCALData: public NCFuncs +{ + std::vector times; + MString datapath; + + struct Parameters: public BaseParameters + { + real lonb, latb, lone, late; + virtual ~Parameters() override = default; + }; + + static MString Title() { return "AVISO local mirror"; } + + public: + using Data = Simple2DData; + + AVISOLOCALData() = default; + + MString Info() const; + // TODO: RetVal + MString Open(const CLArgs& args); + + std::pair Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const; + + bool isOk() const { return times.size() > 0; } + + 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; } + + explicit operator bool() const { return times.size() > 0; } + + bool CheckVar(const MString& vname) const + { + NCFileA nc; + nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc"); + return NCFuncs::CheckVar(vname, [&nc = std::as_const(nc)](const MString& vn) { return HaveVar(nc, vn); }); + } + + Data Read(const MString& vname, const BaseParameters* ip, size_t i) const; + + template Data ReadVarRaw(const NCFileA& nc, const MString& name, const struct Parameters* p) const; + + UVData ReadUV(const BaseParameters* ip, size_t i) const; +}; diff --git a/include/layereddata.h b/include/layereddata.h index 8610fce..e0af873 100644 --- a/include/layereddata.h +++ b/include/layereddata.h @@ -3,6 +3,7 @@ #include "ParseArgs.h" #include "gsw.h" #include "mdatetime.h" +#include "ncfuncs.h" #include "simple2ddata.h" #include "uvdata.h" #include @@ -14,12 +15,13 @@ using michlib::DetGeoDomain; using michlib::Floor; using michlib::GPL; using michlib::int2; +using michlib::int4; using michlib::MDateTime; using michlib::MString; using michlib::NCFileA; using michlib::ToGeoDomain; -class LayeredData +class LayeredData: public NCFuncs { public: using Data = Simple2DData; @@ -34,32 +36,22 @@ class LayeredData public: NC(MString&& newurl): url(std::move(newurl)) { nc.Reset(url + "#cache&noprefetch"); } - MString ReadTimes() + MString ReadTimes(const MString& tname) { if(!nc) return "File not open"; - auto time = nc.VR("time"); + auto time = nc.VR(tname); if(!time) return "Can't read times"; MDateTime refdate; time_t step = 0; { - auto units = nc.Attribute("time", "units"); + auto units = nc.Attribute(tname, "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()) - { - if(*ci == "hours") step = 3600; - if(*ci == "days") step = 3600 * 24; - ci++; - } - 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"; + auto [rd, st, suc] = Refdate(units); + if(!suc) return "Can't parse " + units.Get() + " to refdate"; + if(st == 0) return "Can't get timestep from string " + units.Get(); + refdate = rd; + step = st; } times.resize(time.DimLen(0)); @@ -72,6 +64,8 @@ class LayeredData MDateTime Begin() const { return times.front(); } MDateTime End() const { return times.back(); } + const NCFileA& Get() const { return nc; } + const std::vector& Times() const { return times; } size_t Index(MDateTime tm) const { @@ -94,8 +88,7 @@ class LayeredData std::vector nc; std::vector depths; std::vector times; - MString lonname, latname; - size_t nx, ny; + struct CoordNames dname; real lonb, latb, lone, late; real lonstep, latstep; MString title; @@ -194,16 +187,7 @@ class LayeredData bool CheckVar(const MString& vname) const { - if(!HaveVar(vname)) - { - bool varexist = false; - if(vname == "temp" && HaveVar("ptemp") && HaveVar("sal")) varexist = true; - if(vname == "ptemp" && HaveVar("temp") && HaveVar("sal")) varexist = true; - if(vname == "pdens" && (HaveVar("ptemp") || HaveVar("temp")) && HaveVar("sal")) varexist = true; - if((vname == "U" || vname == "U2") && HaveVar("u") && HaveVar("v")) varexist = true; - if(!varexist) return false; - } - return true; + return NCFuncs::CheckVar(vname, [this](const MString& vn) { return HaveVar(vn); }); } private: @@ -212,15 +196,7 @@ class LayeredData 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; - } - } + if(NCFuncs::HaveVar(nc[i].Get(), vname)) return true; return false; } @@ -240,20 +216,4 @@ class LayeredData } return {"", 0, 0}; } - - static MString StName2Name(const MString& stname) - { - if(stname == "sea_water_potential_temperature") return "ptemp"; - 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 == "eastward_sea_water_velocity") return "u"; - if(stname == "northward_sea_water_velocity") return "v"; - if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "u"; - if(stname == "surface_geostrophic_northward_sea_water_velocity") return "v"; - if(stname == "upward_sea_water_velocity") return "w"; - return ""; - } }; diff --git a/include/ncfuncs.h b/include/ncfuncs.h new file mode 100644 index 0000000..f9445b4 --- /dev/null +++ b/include/ncfuncs.h @@ -0,0 +1,40 @@ +#pragma once +#include "DataAdapters/ncfilealt.h" +#include "StringFunctions.h" +#include "mdatetime.h" +#include +#include + +using michlib::MDateTime; +using michlib::MString; +using michlib::NCFileA; + +class NCFuncs +{ + public: + struct CoordNames + { + MString lonname, latname, depthname, timename; + size_t nx, ny, nz, nt; + }; + + static MString StName2Name(const MString& stname); + static std::tuple Refdate(const MString& refdate); + static void GetVars(const NCFileA& nc, std::set& vars); + static CoordNames GetCNames(const NCFileA& nc); + static CoordNames GetDNames(const NCFileA& nc); + static bool HaveVar(const NCFileA& nc, const MString& vname); + template static bool CheckVar(const MString& vname, HV hv) + { + if(!hv(vname)) + { + bool varexist = false; + if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = true; + if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = true; + if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = true; + if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = true; + if(!varexist) return false; + } + return true; + } +}; diff --git a/include/odm.h b/include/odm.h index 247f9f3..3d73836 100644 --- a/include/odm.h +++ b/include/odm.h @@ -1,5 +1,6 @@ #pragma once #include "AVISO.h" +#include "AVISOLOCAL.h" #include "BFileW.h" #include "HYCOM.h" #include "NEMO.h" @@ -11,7 +12,7 @@ using michlib::errmessage; using michlib::GPL; using michlib::message; -using DataVariants = std::variant; +using DataVariants = std::variant; template concept InfoSupported = requires(T t) { @@ -154,6 +155,13 @@ class Data: public DataVariants if(res.Exist()) return "Can't open source " + src + ":\n" + res; *this = Data(std::move(data)); } + else if(src == "AVISOLOCAL") + { + AVISOLOCALData data; + auto res = data.Open(args); + if(res.Exist()) return "Can't open source " + src + ":\n" + res; + *this = Data(std::move(data)); + } else if(src == "HYCOM") { HYCOMData data; diff --git a/include/simple2ddata.h b/include/simple2ddata.h index 32471bd..2d625ad 100644 --- a/include/simple2ddata.h +++ b/include/simple2ddata.h @@ -43,4 +43,7 @@ class Simple2DData: public BaseData real Lon(size_t i) const { return Lon(i % nx, i / nx); } real Lat(size_t i) const { return Lat(i % nx, i / nx); } + + real XStep() const { return xstep; } + real YStep() const { return ystep; } }; diff --git a/src/AVISOLOCAL.cpp b/src/AVISOLOCAL.cpp new file mode 100644 index 0000000..615682e --- /dev/null +++ b/src/AVISOLOCAL.cpp @@ -0,0 +1,258 @@ +#define MICHLIB_NOSOURCE +#include "AVISOLOCAL.h" + +MString AVISOLOCALData::Info() const +{ + if(!isOk()) return ""; + + NCFileA nc; + struct CoordNames cn, dn; + std::set vars; + + nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc"); + if(!nc) return "Can't open file " + datapath + "/uv-" + times[0].ToString() + ".nc"; + + GetVars(nc, vars); + dn = GetDNames(nc); + cn = GetCNames(nc); + + MString svars; + { + bool first = true; + for(const auto& v: vars) + { + svars += (first ? "" : ", ") + v; + first = false; + } + } + + auto lons = nc.VR(cn.lonname); + auto lats = nc.VR(cn.latname); + if(!(lons && lats)) return "Can't get longitudes/latitudes"; + + real lonb = lons(0); + real latb = lats(0); + real lone = lons(dn.nx - 1); + real late = lats(dn.ny - 1); + real lonstep = (lone - lonb) / (dn.nx - 1); + real latstep = (late - latb) / (dn.ny - 1); + + // 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: " + dn.nx + "x" + dn.ny + " (" + lonstep + " x " + latstep + ")\n" + + " Supported variables: " + svars; + // clang-format on +} + +MString AVISOLOCALData::Open(const CLArgs& args) +{ + GPL.UsePrefix("AVISOLOCAL"); + datapath = GPL.ParameterSValue("Datapath", ""); + RegExp regex("uv-([0-9]{4}-[0-9]{2}-[0-9]{2}).nc"); + regex.Compile(); + + DIR* dir = opendir(datapath.Buf()); + struct dirent* de; + if(nullptr == dir) return "Can't open directory " + datapath; + while((de = readdir(dir))) + { + if(!regex.Match(de->d_name)) continue; + times.emplace_back(MString(de->d_name + regex.Off(1), regex.Len(1))); + } + closedir(dir); + std::sort(times.begin(), times.end()); + return ""; +} + +std::pair AVISOLOCALData::Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const +{ + std::unique_ptr ppar{new struct Parameters}; + if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"}; + ppar->lonb = args.at("lonb").ToReal(); + ppar->lone = args.at("lone").ToReal(); + ppar->latb = args.at("latb").ToReal(); + ppar->late = args.at("late").ToReal(); + + pars.SetParameter("lonb", ppar->lonb); + pars.SetParameter("latb", ppar->latb); + pars.SetParameter("lone", ppar->lone); + pars.SetParameter("late", ppar->late); + + return {ppar.release(), ""}; +} + +AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParameters* ip, size_t i) const +{ + if(!isOk()) return Data(); + + auto p = dynamic_cast(ip); + NCFileA nc; + MString name = ""; + bool isfloat = false, isint2 = false, isint4 = false; + + nc.Reset(datapath + "/uv-" + times[i].ToString() + ".nc"); + if(!nc) return Data(); + { + auto head = nc.Header(); + for(const auto& v: head.Variables()) + { + auto stname = nc.A(v.Name(), "standard_name"); + if(!stname) continue; + if(StName2Name(stname) == vname) + { + name = v.Name(); + isint2 = v.Type().Id() == NC_SHORT; + isint4 = v.Type().Id() == NC_INT; + isfloat = v.Type().Id() == NC_FLOAT; + } + } + } + + if(!name.Exist()) // Conversion read + { + // U and U2 from u and v + if(vname == "U" || vname == "U2") + { + bool square = vname == "U2"; + auto u = Read("u", ip, i); + auto v = Read("v", ip, i); + if(!(u && v)) return Data(); + auto out = u; + for(size_t ind = 0; ind < out.N(); ind++) + { + if(u.IsFill(ind) || v.IsFill(ind)) + out.V(ind) = out.Fillval(); + else + out.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind)); + } + return out; + } + return Data(); + } + + // Direct read + if(isint2) return ReadVarRaw(nc, name, p); + if(isint4) return ReadVarRaw(nc, name, p); + if(isfloat) return ReadVarRaw(nc, name, p); + return Data(); +} + +template AVISOLOCALData::Data AVISOLOCALData::ReadVarRaw(const NCFileA& nc, const MString& name, const struct AVISOLOCALData::Parameters* p) const +{ + real unitmul = 1.0; + DataType fill; + real offset = 0.0, scale = 1.0; + + { + auto a_fill = nc.A(name, "_FillValue"); + auto a_offset_d = nc.A(name, "add_offset"); + auto a_scale_d = nc.A(name, "scale_factor"); + auto a_offset_f = nc.A(name, "add_offset"); + auto a_scale_f = nc.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 cn = GetCNames(nc); + auto dn = GetDNames(nc); + + auto lons = nc.VR(cn.lonname); + auto lats = nc.VR(cn.latname); + if(!(lons && lats)) return Data(); + + auto lonb = lons(0); + auto latb = lats(0); + auto lone = lons(dn.nx - 1); + auto late = lats(dn.ny - 1); + auto lonstep = (lone - lonb) / (dn.nx - 1); + auto latstep = (late - latb) / (dn.ny - 1); + + auto dom = DetGeoDomain(lonb, lone); + real lon1 = ToGeoDomain(p->lonb, dom); + real lon2 = ToGeoDomain(p->lone, dom); + real lat1 = p->latb; + real lat2 = p->late; + + auto yb = static_cast(Floor((lat1 - latb) / latstep)); + auto ye = static_cast(Ceil((lat2 - latb) / latstep)); + if(ye > dn.ny - 1) ye = dn.ny - 1; + if(yb >= ye) return Data(); + + auto xb = static_cast(Floor((lon1 - lonb) / lonstep)); + auto xe = static_cast(Ceil((lon2 - lonb) / lonstep)); + + if(xb == xe) return Data(); + + auto unit = nc.A(name, "units"); + if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0; + + Data data((xb < xe) ? (xe - xb + 1) : (dn.nx + xe - xb + 1), ye - yb + 1, lons(xb), lats(yb), lonstep, latstep); + + if(xb < xe) + { + auto var = nc.V(name, {dn.lonname, xb, xe - xb + 1}, {dn.latname, yb, ye - yb + 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 = nc.V(name, {dn.lonname, xb}, {dn.latname, yb, ye - yb + 1}); + auto var2 = nc.V(name, {dn.lonname, 0, xe + 1}, {dn.latname, yb, ye - yb + 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; +} + +UVData AVISOLOCALData::ReadUV(const BaseParameters* ip, size_t i) const +{ + if(!isOk()) return UVData(); + + auto u = Read("u", ip, i); + auto v = Read("v", ip, i); + if(!(u && v)) return UVData(); + + UVData out{u.Nx(), u.Ny(), u.Lon(0, 0), u.Lat(0, 0), u.XStep(), u.YStep()}; + for(size_t i = 0; i < out.N(); i++) + { + if(u(i) == Data::Fillval() || v(i) == Data::Fillval()) + out.U(i) = out.V(i) = UVData::Fillval(); + else + { + out.U(i) = u(i); + out.V(i) = v(i); + } + } + return out; +} diff --git a/src/layereddata.cpp b/src/layereddata.cpp index 93878a5..991a06c 100644 --- a/src/layereddata.cpp +++ b/src/layereddata.cpp @@ -8,22 +8,7 @@ MString LayeredData::Info() const 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; - if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret)); - } - } - if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens"); - if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp"); - if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp"); - - if(vars.contains("u") && vars.contains("v")) vars.emplace("U"); - if(vars.contains("u") && vars.contains("v")) vars.emplace("U2"); + for(const auto& f: nc) GetVars(f.Get(), vars); MString svars; { @@ -43,7 +28,7 @@ MString LayeredData::Info() const " Time step: " + Timestep() + " seconds\n" + " Time moments: " + NTimes() + "\n" + " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" + - " Grid: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" + + " Grid: " + dname.nx + "x" + dname.ny + " (" + lonstep + " x " + latstep + ")\n" + " Depths:" + d + "\n" + " Supported variables: " + svars; // clang-format on @@ -69,12 +54,6 @@ MString LayeredData::Open(const MString& dataset) 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; @@ -82,60 +61,66 @@ MString LayeredData::Open(const MString& dataset) } if(nc.size() == 0) return "No urls for dataset " + dataset + " specified in config"; - auto head = nc[0]->Header(); - lonname = latname = ""; - for(const auto& dim: head.Dimensions()) + dname = GetDNames(nc[0].Get()); + if(!(dname.lonname.Exist() && dname.latname.Exist())) { - 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(); - } + nc.clear(); + return "Can't find longitude/latitude"; } - if(!(lonname.Exist() && latname.Exist())) + if(!dname.timename.Exist()) { nc.clear(); - return "Can't find longitude/latitude"; + return "Can't find time"; } + auto cn = GetCNames(nc[0].Get()); + // Read times + for(auto& f: nc) { - 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()); + MString ret = f.ReadTimes(cn.timename); + if(ret.Exist()) + { + nc.clear(); + return ret; + } + 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) + if(cn.depthname.Exist()) { + auto rdepths = nc[0]->VR(cn.depthname); + 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); } else // Surface only data { depths.resize(1); - depths[0]=0; + depths[0] = 0; } - auto lons = nc[0]->VR(lonname); - auto lats = nc[0]->VR(latname); + auto lons = nc[0]->VR(cn.lonname); + auto lats = nc[0]->VR(cn.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); + lone = lons(dname.nx - 1); + late = lats(dname.ny - 1); + lonstep = (lone - lonb) / (dname.nx - 1); + latstep = (late - latb) / (dname.ny - 1); return ""; } @@ -144,7 +129,7 @@ std::pair LayeredData::Parameters(michlib_intern { std::unique_ptr ppar{new struct Parameters}; - if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger(); + ppar->layer = args.contains("layer") ? args.at("layer").ToInteger() : 0; 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); @@ -159,7 +144,7 @@ std::pair LayeredData::Parameters(michlib_intern 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->ye > dname.ny - 1) ppar->ye = dname.ny - 1; if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"}; ppar->xb = static_cast(Floor((lon1 - lonb) / lonstep)); @@ -277,7 +262,9 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters* { if(v.Dimensions().size() == 3) nodepth = true; if(v.Type().Id() == NC_SHORT) return ReadVarRaw(nc[id], name, tid, nodepth, p); + if(v.Type().Id() == NC_INT) return ReadVarRaw(nc[id], name, tid, nodepth, p); if(v.Type().Id() == NC_FLOAT) return ReadVarRaw(nc[id], name, tid, nodepth, p); + if(v.Type().Id() == NC_DOUBLE) return ReadVarRaw(nc[id], name, tid, nodepth, p); } return Data(); @@ -317,8 +304,15 @@ template LayeredData::Data LayeredData::ReadVarRaw(const NC& f, 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_fill) + fill = a_fill; + else + { + if constexpr(std::is_floating_point_v) + fill = NAN; + else + fill = -1; + } if(a_offset_d) offset = a_offset_d; if(a_scale_d) scale = a_scale_d; if(a_offset_f) offset = a_offset_f; @@ -328,12 +322,13 @@ template LayeredData::Data LayeredData::ReadVarRaw(const NC& 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); + Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (dname.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, {dname.lonname, p->xb, p->xe - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}) + : f->V(name, {dname.lonname, p->xb, p->xe - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, p->layer, 1}); if(!var) return Data(); if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data(); @@ -341,28 +336,28 @@ template LayeredData::Data LayeredData::ReadVarRaw(const NC& f, 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); + data(ix, iy) = (v == fill || isnan(v)) ? 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}); + auto var1 = nodepth ? f->V(name, {dname.lonname, p->xb}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}) + : f->V(name, {dname.lonname, p->xb}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, p->layer, 1}); + auto var2 = nodepth ? f->V(name, {dname.lonname, 0, p->xe + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}) + : f->V(name, {dname.lonname, 0, p->xe + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, 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); + data(ix, iy) = (v == fill || isnan(v)) ? 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); + data(ix + var1.DimLen(0), iy) = (v == fill || isnan(v)) ? Data::Fillval() : ((v * scale + offset) * unitmul); } } return data; diff --git a/src/ncfuncs.cpp b/src/ncfuncs.cpp new file mode 100644 index 0000000..6110e5c --- /dev/null +++ b/src/ncfuncs.cpp @@ -0,0 +1,148 @@ +#define MICHLIB_NOSOURCE +#include "ncfuncs.h" + +NCFuncs::CoordNames NCFuncs::GetDNames(const NCFileA& nc) +{ + CoordNames out; + auto head = nc.Header(); + for(const auto& dim: head.Dimensions()) + { + if(dim.Name() == "lon" || dim.Name() == "longitude") + { + out.lonname = dim.Name(); + out.nx = dim.Len(); + } + if(dim.Name() == "lat" || dim.Name() == "latitude") + { + out.latname = dim.Name(); + out.ny = dim.Len(); + } + if(dim.Name() == "depth") + { + out.depthname = dim.Name(); + out.nz = dim.Len(); + } + if(dim.Name() == "time") + { + out.timename = dim.Name(); + out.nt = dim.Len(); + } + } + return out; +} + +NCFuncs::CoordNames NCFuncs::GetCNames(const NCFileA& nc) +{ + CoordNames out; + auto head = nc.Header(); + for(const auto& v: head.Variables()) // Try to define coordinates by attribute standard_name or attribute axis + { + auto stname = nc.A(v.Name(), "standard_name"); + auto axis = nc.A(v.Name(), "axis"); + bool islon = false, islat = false, isdepth = false, istime = false; + if(!(stname || axis)) continue; + if(stname && stname.Get() == "longitude") islon = true; + if(stname && stname.Get() == "latitude") islat = true; + if(stname && stname.Get() == "depth") isdepth = true; + if(stname && stname.Get() == "time") istime = true; + + if(!out.lonname.Exist() && axis && axis.Get() == "X") islon = true; + if(!out.latname.Exist() && axis && axis.Get() == "Y") islat = true; + if(!out.depthname.Exist() && axis && axis.Get() == "Z") isdepth = true; + if(!out.timename.Exist() && axis && axis.Get() == "T") istime = true; + + if(islon) out.lonname = v.Name(); + if(islat) out.latname = v.Name(); + if(isdepth) out.depthname = v.Name(); + if(istime) out.timename = v.Name(); + + if(islon) out.nx = v.Dimensions().size(); + if(islat) out.ny = v.Dimensions().size(); + if(isdepth) out.nz = v.Dimensions().size(); + if(istime) out.nt = v.Dimensions().size(); + } + return out; +} + +void NCFuncs::GetVars(const NCFileA& nc, std::set& vars) +{ + auto head = nc.Header(); + for(const auto& v: head.Variables()) + { + auto ret = nc.A(v.Name(), "standard_name"); + if(!ret) continue; + if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret)); + } + if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens"); + if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp"); + if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp"); + + if(vars.contains("u") && vars.contains("v")) vars.emplace("U"); + if(vars.contains("u") && vars.contains("v")) vars.emplace("U2"); +} + +std::tuple NCFuncs::Refdate(const MString& refdate) +{ + MDateTime out; + time_t step = 0; + + MString rstr; + auto words = michlib::Split_on_words(refdate); + auto ci = words.begin(); + if(ci != words.end()) + { + if(*ci == "hours") step = 3600; + if(*ci == "days") step = 3600 * 24; + ci++; + } + 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 + bool success = out.FromString(rstr); + return {out, step, success}; +} + +MString NCFuncs::StName2Name(const MString& stname) +{ + if(stname == "sea_water_potential_temperature") return "ptemp"; + 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 == "eastward_sea_water_velocity") return "u"; + if(stname == "northward_sea_water_velocity") return "v"; + if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "u"; + if(stname == "surface_geostrophic_northward_sea_water_velocity") return "v"; + if(stname == "upward_sea_water_velocity") return "w"; + return ""; +} + +bool NCFuncs::HaveVar(const NCFileA& nc, const MString& vname) +{ + auto head = nc.Header(); + for(const auto& v: head.Variables()) + { + auto stname = nc.A(v.Name(), "standard_name"); + if(!stname) continue; + if(StName2Name(stname) == vname) return true; + } + return false; +} + +/* +template bool NCFuncs::CheckVar(const MString& vname, HV hv) +{ + if(!hv(vname)) + { + bool varexist = false; + if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = true; + if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = true; + if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = true; + if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = true; + if(!varexist) return false; + } + return true; +} +*/ diff --git a/src/odm.cpp b/src/odm.cpp index fddc996..66912ac 100644 --- a/src/odm.cpp +++ b/src/odm.cpp @@ -6,7 +6,7 @@ inline void Usage(const MString& arg0) message("Keys are:"); message(" action. What the program should do. May be: info, tsc, uv. Default: tsc."); message(" Keys for action=info. Print some information about dataset."); - message(" source. Required. May be: NEMO, HYCOM, AVISO"); + message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL"); message(" Keys for source=NEMO"); message(" dataset. Can be DT, NRT or NRT6. Default: DT"); message(" Keys for source=HYCOM"); @@ -14,7 +14,7 @@ inline void Usage(const MString& arg0) message(" Keys for source=AVISO"); message(" dataset. Can be DT, NRT, EckmanDT or EckmanNRT. Default: DT"); message(" Keys for action=tsc. Get temperature, salinity, chlorofill from dataset."); - message(" source. Required. May be: NEMO, HYCOM, AVISO"); + message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL"); message(" var. Required. May be: U, U2, u, v, temp, ptemp, pdens, 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"); @@ -30,7 +30,7 @@ inline void Usage(const MString& arg0) message(" dataset. Can be DT, NRT, EckmanDT or EckmanNRT. Default: DT"); message(" layer and/or depth. Layer or depth of AVISO dataset. If depth is specified, layer is ignored. Both ignored for datasets DT and NRT. Default: layer=0"); message(" Keys for action=uv. Get velocity field and its derivatives."); - message(" source. Required. May be: NEMO, HYCOM, AVISO"); + message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL"); 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 for components of velocity field, divergency, rotor and Okubo-Weiss parameter. If absent, this data not calculated."); @@ -68,10 +68,7 @@ int main(int argc, char** argv) else if(action == "tsc") ret = data.ActionTsc(args); else if(action == "uv") - { - args["var"] = "U"; - ret = data.ActionUV(args); - } + ret = data.ActionUV(args); else { errmessage("Unknown action " + action);