From 04f45ba459d8fd5dd1d0a5ef12b7d8a0b3ef8885 Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Sun, 28 Jan 2024 18:12:06 +1000 Subject: [PATCH] Support for time column in netcdf files generated by uv and tsc actions --- actions/actiontsc.h | 105 ++++++++------ actions/actionuv.cpp | 56 ++++++-- actions/actionuv.h | 323 ++++++++++++++++++++++++------------------- include/ncfilew.h | 33 ++++- src/ncfilew.cpp | 38 ++++- 5 files changed, 350 insertions(+), 205 deletions(-) diff --git a/actions/actiontsc.h b/actions/actiontsc.h index a09312b..f34bdaf 100644 --- a/actions/actiontsc.h +++ b/actions/actiontsc.h @@ -41,6 +41,11 @@ template MString ActionTSC::DoAction(const CLArgs& args, D& ds) } } + bool average = args.contains("average"); + + int compress = 3; + if(args.contains("compress")) compress = args.at("compress").ToInt(); + auto vlist = michlib::Split_on_words(varstring, " ,", false); std::vector vnames{std::make_move_iterator(std::begin(vlist)), std::make_move_iterator(std::end(vlist))}; { @@ -73,60 +78,68 @@ template MString ActionTSC::DoAction(const CLArgs& args, D& ds) } auto p = sourcepars.get(); - auto data = Read(ds, vnames, p, tindexes); - if(data.size() != vnames.size()) return "Can't read data"; - - for(size_t i = 0; i < vnames.size(); i++) - { - if(!data[i].Unit().Exist()) michlib::errmessage("Unknown measurement unit for variable " + vnames[i] + "!"); - if(!data[i].StandartName().Exist()) data[i].SetStandartName(NCFuncs::Name2StName(vnames[i])); - if(!data[i].LongName().Exist()) data[i].SetLongName(NCFuncs::Name2LongName(vnames[i])); - } + TimeData tdata(ds, tindexes); MString name = args.contains("out") ? args.at("out") : "out.bin"; MString outfmt = args.contains("format") ? args.at("format") : (GetExt(name) == "nc" ? "nc" : "bin"); - if(outfmt == "bin") + if(outfmt == "bin" && !average && tindexes.size() > 1) return "Multiple time moments does'nt supported by the bin format"; + + bool headwrited = false; + NCFileW ncfw; + + for(size_t it = 0; it < tindexes.size(); it++) { - BFileW fw; - - size_t ic = 0; - fw.Create(name, 2 + data.size()); - fw.SetColumnName(++ic, "Longitude"); - fw.SetColumnName(++ic, "Latitude"); - for(size_t i = 0; i < data.size(); i++) fw.SetColumnName(++ic, vnames[i] + ", " + (data[i].Unit().Exist() ? data[i].Unit() : "unknown")); - fw.SetParameters(pars); - for(size_t r = 0; r < data[0].N(); r++) + auto data = average ? Read(ds, vnames, p, tindexes) : Read(ds, vnames, p, tindexes[it]); + if(data.size() != vnames.size()) return "Can't read data"; + if(!headwrited) + for(size_t i = 0; i < vnames.size(); i++) + { + if(!data[i].Unit().Exist()) michlib::errmessage("Unknown measurement unit for variable " + vnames[i] + "!"); + if(!data[i].StandartName().Exist()) data[i].SetStandartName(NCFuncs::Name2StName(vnames[i])); + if(!data[i].LongName().Exist()) data[i].SetLongName(NCFuncs::Name2LongName(vnames[i])); + } + + if(outfmt == "bin") { - fw.Write(data[0].Lon(r)); - fw.Write(data[0].Lat(r)); - for(size_t i = 0; i < data.size(); i++) fw.Write(data[i].IsFill(r) ? NAN : data[i](r)); + BFileW fw; + + size_t ic = 0; + fw.Create(name, 2 + data.size()); + fw.SetColumnName(++ic, "Longitude"); + fw.SetColumnName(++ic, "Latitude"); + for(size_t i = 0; i < data.size(); i++) fw.SetColumnName(++ic, vnames[i] + ", " + (data[i].Unit().Exist() ? data[i].Unit() : "unknown")); + fw.SetParameters(pars); + for(size_t r = 0; r < data[0].N(); r++) + { + fw.Write(data[0].Lon(r)); + fw.Write(data[0].Lat(r)); + for(size_t i = 0; i < data.size(); i++) fw.Write(data[i].IsFill(r) ? NAN : data[i](r)); + } + fw.Finalize(); + fw.Close(); } - fw.Finalize(); - fw.Close(); - return ""; - } + else if(outfmt == "nc" || outfmt == "netcdf") + { + MString err; + + if(!err.Exist() && !headwrited) err = ncfw.Create(data[0], name, compress); + if(!err.Exist() && !headwrited) err = ncfw.AddTimeData(tdata, !average); + if(!err.Exist() && !headwrited) err = ncfw.AddAtts(pars); + for(size_t i = 0; i < data.size(); i++) + if(!err.Exist() && !headwrited) err = ncfw.AddVariable(vnames[i], data[i].StandartName(), data[i].LongName(), data[i].Unit(), data[i].Comment()); + if(!err.Exist() && !headwrited) err = ncfw.WriteGrid(data[0]); + for(size_t i = 0; i < data.size(); i++) + if(!err.Exist()) err = average ? ncfw.WriteVariable(data[i], vnames[i]) : ncfw.WriteVariable(data[i], vnames[i], it); + if(err.Exist()) return err; + } + else + return "Unknown format: " + outfmt; - if(outfmt == "nc" || outfmt == "netcdf") - { - int compress = 3; - NCFileW fw; - MString err; - - if(args.contains("compress")) compress = args.at("compress").ToInt(); - - if(!err.Exist()) err = fw.Create(data[0], name, compress); - if(!err.Exist()) err = fw.AddAtts(pars); - for(size_t i = 0; i < data.size(); i++) - if(!err.Exist()) err = fw.AddVariable(vnames[i], data[i].StandartName(), data[i].LongName(), data[i].Unit(), data[i].Comment()); - if(!err.Exist()) err = fw.WriteGrid(data[0]); - for(size_t i = 0; i < data.size(); i++) - if(!err.Exist()) err = fw.WriteVariable(data[i], vnames[i]); - if(err.Exist()) return err; - fw.Close(); - - return ""; + if(average) break; + headwrited = true; } + ncfw.Close(); - return "Unknown format: " + outfmt; + return ""; }; diff --git a/actions/actionuv.cpp b/actions/actionuv.cpp index 92a6261..1c5ffd2 100644 --- a/actions/actionuv.cpp +++ b/actions/actionuv.cpp @@ -20,9 +20,9 @@ void UVMethods::StPoints::WriteBinBile(const MString& name, const michlib_intern stp.Close(); } -MString UVMethods::StPoints::WriteNcFile(const MString& name, const michlib_internal::ParameterListEx& pars, const MString& history, int comp) const +MString UVMethods::StPoints::CreateNcFile(const MString& name, const michlib_internal::ParameterListEx& pars, const MString& history, int comp, const TimeData& tdata, bool timedep) { - NCFileWBase nc; + tdep = timedep; const MString xname = lonlat ? "longitude" : "x"; const MString yname = lonlat ? "latitude" : "y"; @@ -33,8 +33,24 @@ MString UVMethods::StPoints::WriteNcFile(const MString& name, const michlib_inte nc.AddAtt("history", history); nc.AddAtts(pars); nc.AddDim("i", N()); - nc.AddVar(xname, NC_FLOAT, "i"); - nc.AddVar(yname, NC_FLOAT, "i"); + nc.AddDim("time", tdata.steps.size()); + + nc.AddVar("time", decltype(nc)::Type2NCType, "time"); + nc.SetComp("time", comp); + nc.AddAtt("time", "standard_name", "time"); + nc.AddAtt("time", "long_name", "time"); + nc.AddAtt("time", "units", tdata.UnitName()); + + if(tdep) + { + nc.AddVar(xname, NC_FLOAT, "time", "i"); + nc.AddVar(yname, NC_FLOAT, "time", "i"); + } + else + { + nc.AddVar(xname, NC_FLOAT, "i"); + nc.AddVar(yname, NC_FLOAT, "i"); + } nc.SetComp(xname, comp); nc.SetComp(yname, comp); @@ -51,20 +67,38 @@ MString UVMethods::StPoints::WriteNcFile(const MString& name, const michlib_inte nc.AddAtt(yname, "long_name", "y-coordinate"); } - nc.AddVar("type", NC_UBYTE, "i"); + if(tdep) + nc.AddVar("type", NC_UBYTE, "time", "i"); + else + nc.AddVar("type", NC_UBYTE, "i"); nc.SetComp("type", comp); nc.AddAtt("type", "long_name", "Stationary point type, 0 - saddle, 1 - st. anticicl. focus, 2 - st. knot, 3 - unst. anticicl. focus, 4 - unst. knot, 5 - st. cicl. focus, 6 - unst. cicl. focus"); - if(!nc) return "Can't set grid in the netcdf file " + name + ": " + nc.ErrMessage(); + nc.WriteVar("time", tdata.steps.data()); - nc.EndDef(); + if(!nc) return "Can't set grid in the netcdf file " + name + ": " + nc.ErrMessage(); + return ""; +} - nc.WriteVar(xname, x.data()); - nc.WriteVar(yname, y.data()); - nc.WriteVar("type", t.data()); +MString UVMethods::StPoints::WriteNcFile(size_t it) +{ + const MString xname = lonlat ? "longitude" : "x"; + const MString yname = lonlat ? "latitude" : "y"; - if(!nc) return "Can't write data to the netcdf file " + name + ": " + nc.ErrMessage(); + if(tdep) + { + nc.WriteVar(xname, it, x.data()); + nc.WriteVar(yname, it, y.data()); + nc.WriteVar("type", it, t.data()); + } + else + { + nc.WriteVar(xname, x.data()); + nc.WriteVar(yname, y.data()); + nc.WriteVar("type", t.data()); + } + if(!nc) return MString("Can't write data to the netcdf file: ") + nc.ErrMessage(); return ""; } diff --git a/actions/actionuv.h b/actions/actionuv.h index d5d74f7..ea671b4 100644 --- a/actions/actionuv.h +++ b/actions/actionuv.h @@ -70,6 +70,8 @@ class UVMethods std::vector x, y; std::vector t; bool lonlat; + NCFileWBase nc; + bool tdep; size_t N() const { return t.size(); } @@ -86,9 +88,19 @@ class UVMethods } } + void Reset() + { + x.clear(); + y.clear(); + t.clear(); + } + void WriteBinBile(const MString& name, const michlib_internal::ParameterListEx& pars) const; - MString WriteNcFile(const MString& name, const michlib_internal::ParameterListEx& pars, const MString& history, int comp) const; + MString CreateNcFile(const MString& name, const michlib_internal::ParameterListEx& pars, const MString& history, int comp, const TimeData& tdata, bool timedep); + MString WriteNcFile(size_t it); + + void CloseNcFile() { nc.Close(); } }; }; @@ -112,6 +124,11 @@ template MString ActionUV::DoAction(const CLArgs& args, D& ds) auto [tindexes, err] = GetTIndexes(ds, args, pars); if(err.Exist()) return err; + bool average = args.contains("average"); + + int compress = 3; + if(args.contains("compress")) compress = args.at("compress").ToInt(); + std::unique_ptr sourcepars; if constexpr(ParametersSupported) { @@ -130,160 +147,188 @@ template MString ActionUV::DoAction(const CLArgs& args, D& ds) } auto p = sourcepars.get(); - auto data = ReadUV(ds, p, tindexes, mode); - if(!data) return "Can't read data"; + TimeData tdata(ds, tindexes); - // Main file - MString name = args.contains("out") ? args.at("out") : ""; - MString outfmt = args.contains("outformat") ? args.at("outformat") : (GetExt(name) == "nc" ? "nc" : "bin"); + MString name = args.contains("out") ? args.at("out") : ""; + MString outfmt = args.contains("outformat") ? args.at("outformat") : (GetExt(name) == "nc" ? "nc" : "bin"); + MString namevel = args.contains("velout") ? args.at("velout") : ""; + MString outfmtvel = args.contains("veloutformat") ? args.at("veloutformat") : (GetExt(name) == "nc" ? "nc" : "bin"); + MString namestp = args.contains("stpout") ? args.at("stpout") : ""; + MString outfmtstp = args.contains("stpoutformat") ? args.at("stpoutformat") : (GetExt(name) == "nc" ? "nc" : "bin"); - MString u = data.Unit().Exist() ? data.Unit() : "unknown", d = data.DUnit().Exist() ? data.DUnit() : "unknown"; + size_t shiftx = args.contains("shiftx") ? args.at("shiftx").ToInteger() : 0; + size_t shifty = args.contains("shifty") ? args.at("shifty").ToInteger() : 0; + size_t skipx = args.contains("skipx") ? args.at("skipx").ToInteger() : 1; + size_t skipy = args.contains("skipy") ? args.at("skipy").ToInteger() : 1; - if(name.Exist()) - { - if(outfmt == "bin") - { - BFileW fw; - fw.Create(name, 9); - fw.SetColumnName(1, "Longitude"); - fw.SetColumnName(2, "Latitude"); - fw.SetColumnName(3, "u, " + u); - fw.SetColumnName(4, "v, " + u); - fw.SetColumnName(5, "Div, (" + u + ")/" + d); - fw.SetColumnName(6, "Rot, (" + u + ")/" + d); - fw.SetColumnName(7, "Okubo-Weiss parameter, (" + u + ")2/" + d + "2"); - fw.SetColumnName(8, "Kinetic energy, (" + u + ")2"); - fw.SetColumnName(9, "Eddy kinetic energy, (" + u + ")2"); - fw.SetParameters(pars); - - for(size_t i = 0; i < data.N(); i++) - { - fw.Write(data.Lon(i)); - fw.Write(data.Lat(i)); - fw.Write(data.U(i) == data.Fillval() ? NAN : data.U(i)); - fw.Write(data.V(i) == data.Fillval() ? NAN : data.V(i)); - fw.Write(data.Div(i) == data.Fillval() ? NAN : data.Div(i)); - fw.Write(data.Rot(i) == data.Fillval() ? NAN : data.Rot(i)); - fw.Write(data.OW(i) == data.Fillval() ? NAN : data.OW(i)); - fw.Write(data.U2(i) == data.Fillval() ? NAN : data.U2(i)); - fw.Write((data.U(i) == data.Fillval() || data.V(i) == data.Fillval()) ? NAN : (data.U2(i) - (data.U(i) * data.U(i) + data.V(i) * data.V(i)))); - } - fw.Finalize(); - fw.Close(); - } - else if(outfmt == "nc" || outfmt == "netcdf") - { - int compress = 3; - NCFileW fw; - MString err; - - if(args.contains("compress")) compress = args.at("compress").ToInt(); - - if(!err.Exist()) err = fw.Create(data, name, compress); - if(!err.Exist()) err = fw.AddAtts(pars); - if(!err.Exist()) err = fw.AddVariable("u", "", "Eastward velocity", u, ""); - if(!err.Exist()) err = fw.AddVariable("v", "", "Northward velocity", u, ""); - if(!err.Exist()) err = fw.AddVariable("div", "", "Velocity divergence", "(" + u + ")/" + d, ""); - if(!err.Exist()) err = fw.AddVariable("rot", "", "Velocity rotor", "(" + u + ")/" + d, ""); - if(!err.Exist()) err = fw.AddVariable("ow", "", "Okubo-Weiss parameter", "(" + u + ")2/" + d + "2", ""); - if(!err.Exist()) err = fw.AddVariable("ke", "", "Squared velocity module, u^2+v^2", "(" + u + ")2", ""); - if(!err.Exist()) err = fw.AddVariable("eke", "", "Squared velocity dispersion aka eddy kinetic energy, -^2-^2", "(" + u + ")2", ""); - - if(!err.Exist()) err = fw.WriteGrid(data); - - if(!err.Exist()) err = fw.WriteVariable(data, "u", [&data = std::as_const(data)](size_t i, size_t j) { return data.U(i, j); }); - if(!err.Exist()) err = fw.WriteVariable(data, "v", [&data = std::as_const(data)](size_t i, size_t j) { return data.V(i, j); }); - if(!err.Exist()) err = fw.WriteVariable(data, "div", [&data = std::as_const(data)](size_t i, size_t j) { return data.Div(i, j); }); - if(!err.Exist()) err = fw.WriteVariable(data, "rot", [&data = std::as_const(data)](size_t i, size_t j) { return data.Rot(i, j); }); - if(!err.Exist()) err = fw.WriteVariable(data, "ow", [&data = std::as_const(data)](size_t i, size_t j) { return data.OW(i, j); }); - if(!err.Exist()) err = fw.WriteVariable(data, "ke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j); }); - if(!err.Exist()) - err = fw.WriteVariable(data, "eke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j) - (data.U(i, j) * data.U(i, j) + data.V(i, j) * data.V(i, j)); }); - if(err.Exist()) return err; + if(!average && tindexes.size() > 1 && ((name.Exist() && outfmt == "bin") || (namevel.Exist() && outfmtvel == "bin") || (namestp.Exist() && outfmtstp == "bin"))) + return "Multiple time moments does'nt supported by the bin format"; - fw.Close(); - } - else - return "Unknown format: " + outfmt; - } + bool headwrited = false; + + MString velunit, distunit; + NCFileW fw, fwfilt; + + StPoints stp(true); - // Filtered vectors file - name = args.contains("velout") ? args.at("velout") : ""; - outfmt = args.contains("veloutformat") ? args.at("veloutformat") : (GetExt(name) == "nc" ? "nc" : "bin"); - if(name.Exist()) + for(size_t it = 0; it < tindexes.size(); it++) { - size_t shiftx = args.contains("shiftx") ? args.at("shiftx").ToInteger() : 0; - size_t shifty = args.contains("shifty") ? args.at("shifty").ToInteger() : 0; - size_t skipx = args.contains("skipx") ? args.at("skipx").ToInteger() : 1; - size_t skipy = args.contains("skipy") ? args.at("skipy").ToInteger() : 1; + auto data = average ? ReadUV(ds, p, tindexes, mode) : ReadUV(ds, p, tindexes[it], mode); + if(!headwrited) + { + velunit = data.Unit().Exist() ? data.Unit() : "unknown"; + distunit = data.DUnit().Exist() ? data.DUnit() : "unknown"; + } - Sparser sdata(data, shiftx, shifty, skipx, skipy); - if(outfmt == "bin") + // Main file + if(name.Exist()) { - BFileW vel; - vel.Create(name, 4); - vel.SetColumnName(1, "Longitude"); - vel.SetColumnName(2, "Latitude"); - vel.SetColumnName(3, "u, " + u); - vel.SetColumnName(4, "v, " + u); - vel.SetParameters(pars); - - for(size_t ix = 0; ix < sdata.Nx(); ix++) - for(size_t iy = 0; iy < sdata.Ny(); iy++) + if(outfmt == "bin") + { + BFileW fw; + fw.Create(name, 9); + fw.SetColumnName(1, "Longitude"); + fw.SetColumnName(2, "Latitude"); + fw.SetColumnName(3, "u, " + velunit); + fw.SetColumnName(4, "v, " + velunit); + fw.SetColumnName(5, "Div, (" + velunit + ")/" + distunit); + fw.SetColumnName(6, "Rot, (" + velunit + ")/" + distunit); + fw.SetColumnName(7, "Okubo-Weiss parameter, (" + velunit + ")2/" + distunit + "2"); + fw.SetColumnName(8, "Kinetic energy, (" + velunit + ")2"); + fw.SetColumnName(9, "Eddy kinetic energy, (" + velunit + ")2"); + fw.SetParameters(pars); + + for(size_t i = 0; i < data.N(); i++) { - vel.Write(sdata.Lon(ix, iy)); - vel.Write(sdata.Lat(ix, iy)); - vel.Write(sdata.U(ix, iy) == sdata.Fillval() ? NAN : sdata.U(ix, iy)); - vel.Write(sdata.V(ix, iy) == sdata.Fillval() ? NAN : sdata.V(ix, iy)); + fw.Write(data.Lon(i)); + fw.Write(data.Lat(i)); + fw.Write(data.U(i) == data.Fillval() ? NAN : data.U(i)); + fw.Write(data.V(i) == data.Fillval() ? NAN : data.V(i)); + fw.Write(data.Div(i) == data.Fillval() ? NAN : data.Div(i)); + fw.Write(data.Rot(i) == data.Fillval() ? NAN : data.Rot(i)); + fw.Write(data.OW(i) == data.Fillval() ? NAN : data.OW(i)); + fw.Write(data.U2(i) == data.Fillval() ? NAN : data.U2(i)); + fw.Write((data.U(i) == data.Fillval() || data.V(i) == data.Fillval()) ? NAN : (data.U2(i) - (data.U(i) * data.U(i) + data.V(i) * data.V(i)))); } - vel.Finalize(); - vel.Close(); + fw.Finalize(); + fw.Close(); + } + else if(outfmt == "nc" || outfmt == "netcdf") + { + MString err; + + if(!err.Exist() && !headwrited) err = fw.Create(data, name, compress); + if(!err.Exist() && !headwrited) err = fw.AddTimeData(tdata, !average); + if(!err.Exist() && !headwrited) err = fw.AddAtts(pars); + if(!err.Exist() && !headwrited) err = fw.AddVariable("u", "", "Eastward velocity", velunit, ""); + if(!err.Exist() && !headwrited) err = fw.AddVariable("v", "", "Northward velocity", velunit, ""); + if(!err.Exist() && !headwrited) err = fw.AddVariable("div", "", "Velocity divergence", "(" + velunit + ")/" + distunit, ""); + if(!err.Exist() && !headwrited) err = fw.AddVariable("rot", "", "Velocity rotor", "(" + velunit + ")/" + distunit, ""); + if(!err.Exist() && !headwrited) err = fw.AddVariable("ow", "", "Okubo-Weiss parameter", "(" + velunit + ")2/" + distunit + "2", ""); + if(!err.Exist() && !headwrited) err = fw.AddVariable("ke", "", "Squared velocity module, u^2+v^2", "(" + velunit + ")2", ""); + if(!err.Exist() && !headwrited) err = fw.AddVariable("eke", "", "Squared velocity dispersion aka eddy kinetic energy, -^2-^2", "(" + velunit + ")2", ""); + if(!err.Exist() && !headwrited) err = fw.WriteGrid(data); + + if(!err.Exist()) + err = fw.WriteVariable( + data, "u", [&data = std::as_const(data)](size_t i, size_t j) { return data.U(i, j); }, it); + if(!err.Exist()) + err = fw.WriteVariable( + data, "v", [&data = std::as_const(data)](size_t i, size_t j) { return data.V(i, j); }, it); + if(!err.Exist()) + err = fw.WriteVariable( + data, "div", [&data = std::as_const(data)](size_t i, size_t j) { return data.Div(i, j); }, it); + if(!err.Exist()) + err = fw.WriteVariable( + data, "rot", [&data = std::as_const(data)](size_t i, size_t j) { return data.Rot(i, j); }, it); + if(!err.Exist()) + err = fw.WriteVariable( + data, "ow", [&data = std::as_const(data)](size_t i, size_t j) { return data.OW(i, j); }, it); + if(!err.Exist()) + err = fw.WriteVariable( + data, "ke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j); }, it); + if(!err.Exist()) + err = fw.WriteVariable( + data, "eke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j) - (data.U(i, j) * data.U(i, j) + data.V(i, j) * data.V(i, j)); }, it); + if(err.Exist()) return err; + } + else + return "Unknown format: " + outfmt; } - else if(outfmt == "nc" || outfmt == "netcdf") - { - int compress = 3; - NCFileW fw; - MString err; - - if(args.contains("compress")) compress = args.at("compress").ToInt(); - - if(!err.Exist()) err = fw.Create(sdata, name, compress); - if(!err.Exist()) err = fw.AddAtts(pars); - if(!err.Exist()) err = fw.AddVariable("u", "", "Eastward velocity", u, ""); - if(!err.Exist()) err = fw.AddVariable("v", "", "Northward velocity", u, ""); - - if(!err.Exist()) err = fw.WriteGrid(sdata); - - if(!err.Exist()) err = fw.WriteVariable(sdata, "u", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.U(i, j); }); - if(!err.Exist()) err = fw.WriteVariable(sdata, "v", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.V(i, j); }); - if(err.Exist()) return err; - fw.Close(); + // Filtered file + if(namevel.Exist()) + { + Sparser sdata(data, shiftx, shifty, skipx, skipy); + if(outfmtvel == "bin") + { + BFileW vel; + vel.Create(name, 4); + vel.SetColumnName(1, "Longitude"); + vel.SetColumnName(2, "Latitude"); + vel.SetColumnName(3, "u, " + velunit); + vel.SetColumnName(4, "v, " + velunit); + vel.SetParameters(pars); + + for(size_t ix = 0; ix < sdata.Nx(); ix++) + for(size_t iy = 0; iy < sdata.Ny(); iy++) + { + vel.Write(sdata.Lon(ix, iy)); + vel.Write(sdata.Lat(ix, iy)); + vel.Write(sdata.U(ix, iy) == sdata.Fillval() ? NAN : sdata.U(ix, iy)); + vel.Write(sdata.V(ix, iy) == sdata.Fillval() ? NAN : sdata.V(ix, iy)); + } + vel.Finalize(); + vel.Close(); + } + else if(outfmtvel == "nc" || outfmtvel == "netcdf") + { + MString err; + + if(!err.Exist() && !headwrited) err = fwfilt.Create(sdata, name, compress); + if(!err.Exist() && !headwrited) err = fwfilt.AddTimeData(tdata, !average); + if(!err.Exist() && !headwrited) err = fwfilt.AddAtts(pars); + if(!err.Exist() && !headwrited) err = fwfilt.AddVariable("u", "", "Eastward velocity", velunit, ""); + if(!err.Exist() && !headwrited) err = fwfilt.AddVariable("v", "", "Northward velocity", velunit, ""); + if(!err.Exist() && !headwrited) err = fwfilt.WriteGrid(sdata); + + if(!err.Exist()) + err = fwfilt.WriteVariable( + sdata, "u", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.U(i, j); }, it); + if(!err.Exist()) + err = fwfilt.WriteVariable( + sdata, "v", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.V(i, j); }, it); + if(err.Exist()) return err; + } + else + return "Unknown format: " + outfmtvel; } - else - return "Unknown format: " + outfmt; - } - - // Stationary points - name = args.contains("stpout") ? args.at("stpout") : ""; - outfmt = args.contains("stpoutformat") ? args.at("stpoutformat") : (GetExt(name) == "nc" ? "nc" : "bin"); - if(name.Exist()) - { - StPoints p(true); - for(size_t ix = 0; ix < data.Nx() - 1; ix++) - for(size_t iy = 0; iy < data.Ny() - 1; iy++) p.Add(data.StablePoints(ix, iy)); - if(outfmt == "bin") - p.WriteBinBile(name, pars); - else if(outfmt == "nc" || outfmt == "netcdf") + // Stationary points + if(namestp.Exist()) { - int compress = args.contains("compress") ? args.at("compress").ToInt() : 3; - MString err = p.WriteNcFile(name, pars, args.at("_cmdline"), compress); - if(err.Exist()) return err; + stp.Reset(); + for(size_t ix = 0; ix < data.Nx() - 1; ix++) + for(size_t iy = 0; iy < data.Ny() - 1; iy++) stp.Add(data.StablePoints(ix, iy)); + + if(outfmtstp == "bin") + stp.WriteBinBile(name, pars); + else if(outfmtstp == "nc" || outfmtstp == "netcdf") + { + MString err; + if(!err.Exist() && !headwrited) err = stp.CreateNcFile(name, pars, args.at("_cmdline"), compress, tdata, !average); + if(!err.Exist()) err = stp.WriteNcFile(it); + if(err.Exist()) return err; + } + else + return "Unknown format: " + outfmt; } - else - return "Unknown format: " + outfmt; - } + if(average) break; + headwrited = true; + } + fw.Close(); + fwfilt.Close(); + stp.CloseNcFile(); return ""; }; diff --git a/include/ncfilew.h b/include/ncfilew.h index bcd5d59..02594b8 100644 --- a/include/ncfilew.h +++ b/include/ncfilew.h @@ -1,5 +1,6 @@ #pragma once #include "MString.h" +#include "actiondep.h" #include "traits.h" #include #include @@ -366,6 +367,7 @@ class NCFileW: public NCFileWBase Type type = UNKNOWN; int compress; + bool tdep = false; template static constexpr Type DetType() { @@ -403,9 +405,12 @@ class NCFileW: public NCFileWBase type = UNKNOWN; } + MString AddTimeData(const TimeData& tdata, bool tisindex); MString AddVariable(const MString& name, const MString& stname, const MString& lname, const MString& units, const MString& comment); - template MString WriteVariable(const D& data, const MString& name, Op op) + template + requires(std::is_invocable_v || std::is_invocable_v) + MString WriteVariable(const D& data, const MString& name, Op op, size_t tind) { static constexpr auto dtype = DetType(); if(type == UNKNOWN) return "File not open"; @@ -419,7 +424,10 @@ class NCFileW: public NCFileWBase const size_t c = data.N(); float buf[c]; for(size_t ix = 0; ix < c; ix++) buf[ix] = data.IsFill(ix) ? fill : op(ix); - WriteVar(name, buf); + if(tdep) + WriteVar(name, tind, buf); + else + WriteVar(name, buf); } else { @@ -427,7 +435,10 @@ class NCFileW: public NCFileWBase float buf[c[0] * c[1]]; for(size_t iy = 0; iy < c[0]; iy++) for(size_t ix = 0; ix < c[1]; ix++) buf[iy * c[1] + ix] = data.IsFill(ix, iy) ? fill : op(ix, iy); - WriteVar(name, buf); + if(tdep) + WriteVar(name, tind, buf); + else + WriteVar(name, buf); } if(!*this) return "Can't write variable " + name + ": " + ErrMessage(); @@ -435,14 +446,24 @@ class NCFileW: public NCFileWBase return ""; } - template MString WriteVariable(const D& data, const MString& name) + template MString WriteVariable(const D& data, const MString& name, size_t tind) { if constexpr(Is1DType(DetType())) - return WriteVariable(data, name, [&data = std::as_const(data)](size_t i) { return data(i); }); + return WriteVariable( + data, name, [&data = std::as_const(data)](size_t i) { return data(i); }, tind); else - return WriteVariable(data, name, [&data = std::as_const(data)](size_t i, size_t j) { return data(i, j); }); + return WriteVariable( + data, name, [&data = std::as_const(data)](size_t i, size_t j) { return data(i, j); }, tind); } + template + requires(std::is_invocable_v || std::is_invocable_v) + MString WriteVariable(const D& data, const MString& name, Op op) + { + return WriteVariable(data, name, op, 0); + } + template MString WriteVariable(const D& data, const MString& name) { return WriteVariable(data, name, 0); } + template MString WriteGrid(const D& data) { static constexpr auto dtype = DetType(); diff --git a/src/ncfilew.cpp b/src/ncfilew.cpp index 07ee85c..4d4c26a 100644 --- a/src/ncfilew.cpp +++ b/src/ncfilew.cpp @@ -87,16 +87,48 @@ MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, int compre return ""; } +MString NCFileW::AddTimeData(const TimeData& tdata, bool tisindex) +{ + tdep = tisindex; + + AddDim("time", tdata.steps.size()); + AddVar("time", Type2NCType, "time"); + SetComp("time", compress); + AddAtt("time", "standard_name", "time"); + AddAtt("time", "long_name", "time"); + AddAtt("time", "units", tdata.UnitName()); + WriteVar("time", tdata.steps.data()); + + if(!*this) return MString("Can't add time data to the netcdf file: ") + ErrMessage(); + + return ""; +} + MString NCFileW::AddVariable(const MString& name, const MString& stname, const MString& lname, const MString& units, const MString& comment) { if(type == UNKNOWN) return "File not opened"; if(Is1DType(type)) - AddVar(name, NC_FLOAT, "i"); + { + if(tdep) + AddVar(name, NC_FLOAT, "time", "i"); + else + AddVar(name, NC_FLOAT, "i"); + } else if(IsGeoType(type)) - AddVar(name, NC_FLOAT, "latitude", "longitude"); + { + if(tdep) + AddVar(name, NC_FLOAT, "time", "latitude", "longitude"); + else + AddVar(name, NC_FLOAT, "latitude", "longitude"); + } else - AddVar(name, NC_FLOAT, "y", "x"); + { + if(tdep) + AddVar(name, NC_FLOAT, "time", "y", "x"); + else + AddVar(name, NC_FLOAT, "y", "x"); + } SetComp(name, compress); if(stname.Exist()) AddAtt(name, "standard_name", stname);