From 8ca3fa7bb26b1e7accd8a5714f628034a9810af4 Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Tue, 16 May 2023 16:06:20 +1000 Subject: [PATCH] Added new action GENINTFILE --- include/actiongenintfile.h | 101 +++++++++++++++++++++++++++++++++++++ include/actions.h | 50 +++++++++++------- include/odm.h | 1 + include/traits.h | 30 ++++++++--- include/uvdata.h | 17 ++----- src/odm.cpp | 8 ++- 6 files changed, 167 insertions(+), 40 deletions(-) create mode 100644 include/actiongenintfile.h diff --git a/include/actiongenintfile.h b/include/actiongenintfile.h new file mode 100644 index 0000000..fa74910 --- /dev/null +++ b/include/actiongenintfile.h @@ -0,0 +1,101 @@ +#pragma once +#include "actions.h" + +using michlib::Cos; + +template struct DoAction_ +{ + static MString DoAction(const CLArgs& args, D& data); +}; + +template MString DoAction_::DoAction(const CLArgs& args, D& ds) +{ + michlib_internal::ParameterListEx pars; + pars.UsePrefix(""); + pars.SetParameter("source", args.at("source")); + + if(!args.contains("out")) return "Output file name not specified"; + + auto [tindexes, err] = GetTIndexes(ds, args, pars); + if(err.Exist()) return err; + if(tindexes.size() < 10) return "Too few time moments"; + // Check isotropy + for(size_t it = 0; it < tindexes.size() - 1; it++) + if(ds.Time(tindexes[it + 1]) - ds.Time(tindexes[it]) != ds.Time(tindexes[1]) - ds.Time(tindexes[0])) return "The time step is not isotropic"; + + std::unique_ptr sourcepars; + if constexpr(ParametersSupported) + { + auto [p, err] = ds.Parameters(pars, args); + if(err.Exist()) return err; + sourcepars.reset(p); + } + auto p = sourcepars.get(); + + MString name = args.at("out"); + BFileW fw; + fw.Create(name, 2); + fw.SetColumnName(1, "u"); + fw.SetColumnName(2, "v"); + + for(size_t it = 0; it < tindexes.size(); it++) + { + auto data = ReadUV(ds, p, tindexes[it]); + if(!data) return "Can't read data"; + if(it == 0) + { + real lonb = data.Lon(0, 0), latb = data.Lat(0, 0), lone = data.Lon(data.Nx() - 1, data.Ny() - 1), late = data.Lat(data.Nx() - 1, data.Ny() - 1); + fw.UsePrefix(""); + fw.SetParameter("nx", data.Nx()); + fw.SetParameter("ny", data.Ny()); + fw.SetParameter("dx", data.XStep() * 60.0); + fw.SetParameter("dy", data.YStep() * 60.0); + fw.SetParameter("nt", tindexes.size()); + fw.SetParameter("dt", (ds.Time(tindexes[1]) - ds.Time(tindexes[0])) / 86400.0); + fw.SetParameter("FillValue", data.Fillval()); + fw.UsePrefix("Info"); + fw.SetParameter("lonb", lonb); + fw.SetParameter("latb", latb); + fw.SetParameter("lone", lone); + fw.SetParameter("late", late); + fw.SetParameter("DataSource", ds.DataTitle()); + fw.SetParameter("xy2lon", MString(lonb) + "+%x/60"); + fw.SetParameter("xy2lon_rp", MString("POP 60 DIV ") + lonb + " ADD"); + fw.SetParameter("lonlat2x", MString("(%lon-") + lonb + ")*60"); + fw.SetParameter("xy2lat", MString(latb) + "+%y/60"); + fw.SetParameter("xy2lat_rp", MString("EXCH POP 60 DIV ") + latb + " ADD"); + fw.SetParameter("lonlat2y", MString("(%lat-") + latb + ")*60"); + fw.SetParameter("lonlatuv2dotx", "%u*(0.864/1.852)/Cos(%lat*M_PI/180)"); + fw.SetParameter("lonlatuv2doty", "%v*(0.864/1.852)"); + fw.SetParameter("dotxdotylonlat2u_rp", "EXCH POP EXCH POP COSD 1.852 0.864 DIV MUL MUL"); + fw.SetParameter("dotxdotylonlat2v_rp", "POP POP 1.852 0.864 DIV MUL EXCH POP"); + fw.SetParameter("dotxdotylonlat2u", "(1.852/0.864)*Cos(%lat*M_PI/180)*%dotx"); + fw.SetParameter("dotxdotylonlat2v", "(1.852/0.864)*%doty"); + fw.SetParameter("BeginDate", ds.Time(tindexes.front()).ToTString()); + fw.SetParameter("EndDate", ds.Time(tindexes.back()).ToTString()); + fw.SetParameter("Timestep", ds.Time(tindexes[1]) - ds.Time(tindexes[0])); + fw.SetParameter("DataID", michlib::UniqueId()); + } + + for(size_t ilat = 0; ilat < data.Ny(); ilat++) + { + for(size_t ilon = 0; ilon < data.Nx(); ilon++) + { + if(data.IsCoast(ilon, ilat)) + { + fw.Write(data.Fillval()); + fw.Write(data.Fillval()); + } + else + { + fw.Write(data.U(ilon, ilat) * (0.864 / 1.852) / Cos(data.Lat(ilon, ilat) * M_PI / 180.0)); + fw.Write(data.V(ilon, ilat) * (0.864 / 1.852)); + } + } + } + } + fw.Finalize(); + fw.Close(); + + return ""; +}; diff --git a/include/actions.h b/include/actions.h index e4e952e..fba0b01 100644 --- a/include/actions.h +++ b/include/actions.h @@ -14,7 +14,8 @@ enum class ActionID { INFO, TSC, - UV + UV, + GENINTFILE }; template struct ActionCarrier @@ -39,10 +40,14 @@ template struct IsActionSupportedTest { static constexpr bool ans = ReadPSupported || ReadSupported; }; +template struct IsActionSupportedTest +{ + static constexpr bool ans = ReadIsGrid; +}; template constexpr bool IsActionSupported = IsActionSupportedTest::ans; -using ActionVariants = std::variant, ActionCarrier, ActionCarrier>; +using ActionVariants = std::variant, ActionCarrier, ActionCarrier, ActionCarrier>; class Action: public ActionVariants { @@ -58,6 +63,8 @@ class Action: public ActionVariants *this = ActionVariants(ActionCarrier()); else if(act == "uv") *this = ActionVariants(ActionCarrier()); + else if(act == "genintfile") + *this = ActionVariants(ActionCarrier()); else return "Unknown action: " + act; return ""; @@ -192,32 +199,37 @@ template BaseData Read(const D& data, const MString& vname, const BaseP return BaseData(); } -template UVData> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex) +template UVData> ReadUV(const D& data, const BaseParameters* p, size_t ind) { using RT = ReadType; using UV = UVData; - size_t ind; - if(tindex.size() == 1) + michlib::message("Time: " + data.Time(ind).ToTString()); + RT u, v; + if constexpr(ReadPSupported) { - ind = tindex[0]; - michlib::message("Time: " + data.Time(ind).ToTString()); - RT u, v; - if constexpr(ReadPSupported) - { - u = data.Read("u", p, ind); - v = data.Read("v", p, ind); - } - else if constexpr(ReadSupported) - { - u = data.Read("u", ind); - v = data.Read("v", ind); - } - return UV(u, v); + u = data.Read("u", p, ind); + v = data.Read("v", p, ind); + } + else if constexpr(ReadSupported) + { + u = data.Read("u", ind); + v = data.Read("v", ind); } + return UV(u, v); +} + +template UVData> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex) +{ + using RT = ReadType; + using UV = UVData; + + if(tindex.size() == 1) + return ReadUV(data, p, tindex[0]); else { Averager out; bool ok = true; + size_t ind; for(size_t i = 0; i < tindex.size(); i++) { if(!ok) break; diff --git a/include/odm.h b/include/odm.h index d9ce488..2a3d75b 100644 --- a/include/odm.h +++ b/include/odm.h @@ -1,4 +1,5 @@ #pragma once +#include "actiongenintfile.h" #include "actioninfo.h" #include "actiontsc.h" #include "actionuv.h" diff --git a/include/traits.h b/include/traits.h index 84bc5b5..0345059 100644 --- a/include/traits.h +++ b/include/traits.h @@ -32,13 +32,6 @@ concept ReadSupported = requires(T t, const MString& vname, size_t i) { } -> std::convertible_to; }; -template -concept ReadUVPSupported = requires(T t, const BaseParameters* ip, size_t i) { - { - t.ReadUV(ip, i).U(0, 0) + t.template ReadUV(ip, i).V(0, 0) - } -> std::convertible_to; -}; - template concept IsUVData = requires(T t) { { @@ -52,6 +45,19 @@ concept IsUVData = requires(T t) { } -> std::convertible_to; }; +namespace internal +{ +template +concept HaveXYStep = requires(T t) { + { + t.XStep() + } -> std::convertible_to; + { + t.YStep() + } -> std::convertible_to; +}; +} + namespace internal { template using void_ = void; @@ -72,3 +78,13 @@ template struct GetReadType_ } // namespace internal template using ReadType = internal::GetReadType_, ReadSupported>::type; + +template +concept ReadIsGrid = requires { + { + ReadType().XStep() + } -> std::convertible_to; + { + ReadType().YStep() + } -> std::convertible_to; +}; diff --git a/include/uvdata.h b/include/uvdata.h index ba358c7..889b066 100644 --- a/include/uvdata.h +++ b/include/uvdata.h @@ -1,24 +1,12 @@ #pragma once #include "geohelpers.h" +#include "traits.h" #include using michlib::M_PI; using michlib::real; using michlib::Sqrt; -namespace internal -{ -template -concept HaveXYStep = requires(T t) { - { - t.XStep() - } -> std::convertible_to; - { - t.YStep() - } -> std::convertible_to; -}; -} - enum class Metric { EUCLID, @@ -147,6 +135,9 @@ template class UVDataDims: public UVDataStor real Lat([[maybe_unused]] size_t ix, size_t iy) const { return y0 + iy * ystep; } 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; } }; template class UVDataDims: public UVDataStorage diff --git a/src/odm.cpp b/src/odm.cpp index 7fbfaa6..149552b 100644 --- a/src/odm.cpp +++ b/src/odm.cpp @@ -4,7 +4,7 @@ inline void Usage(const MString& arg0) { message(arg0 + " (=...)"); message("Keys are:"); - message(" action. What the program should do. May be: info, tsc, uv. Default: info."); + message(" action. What the program should do. May be: info, tsc, uv, genintfile. Default: info."); message(" Keys for action=info. Print some information about dataset."); message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL, BINFILE"); message(" Keys for source=NEMO"); @@ -41,6 +41,12 @@ inline void Usage(const MString& arg0) message(" skipx, skipy. How many poinst skipped in the sparsed grid. Used only if velout parameter is present. Default: 1, output each point"); message(" stpout. Output file for stationary points. If absent, this data not calculated."); message(" Keys for specific sources see in the section action=tsc."); + message(" Keys for action=genintfile. Create file with velocity field for interpolation."); + 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, required."); + message(" Keys for specific sources see in the section action=tsc."); } int main(int argc, char** argv)