From a04556729064b96570cd82ead66c00946c9c0715 Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Fri, 14 Apr 2023 14:12:26 +1000 Subject: [PATCH] Added action uv --- include/basedata.h | 54 ++++++++-- include/layereddata.h | 5 + include/odm.h | 8 +- include/uvdata.h | 225 ++++++++++++++++++++++++++++++++++++++++++ src/actiontsc.cpp | 3 + src/actionuv.cpp | 169 +++++++++++++++++++++++++++++++ src/layereddata.cpp | 22 +++++ src/odm.cpp | 20 +++- 8 files changed, 490 insertions(+), 16 deletions(-) create mode 100644 include/uvdata.h create mode 100644 src/actionuv.cpp diff --git a/include/basedata.h b/include/basedata.h index 3bfac64..2a4e445 100644 --- a/include/basedata.h +++ b/include/basedata.h @@ -49,9 +49,34 @@ class BaseData explicit operator bool() const { return N() != 0; } }; +template +concept HaveU = requires(T t) { + { + t.template U(0) + } -> std::convertible_to; +}; + +template +concept HaveV = requires(T t) { + { + t.template V(0) + } -> std::convertible_to; +}; + template class Averager: public Data { - std::vector count; + std::vector count; + static constexpr bool haveu = HaveU; + static constexpr bool havev = HaveV; + + static bool DataOk(const Data* d, size_t i) + { + bool uok = true; + bool vok = true; + if constexpr(haveu) uok = d->U(i) != Data::Fillval(); + if constexpr(havev) uok = d->V(i) != Data::Fillval(); + return uok && vok; + } public: void Add(const Data& d) @@ -61,23 +86,25 @@ template class Averager: public Data { *static_cast(this) = d; count.resize(Data::N()); - for(size_t i = 0; i < Data::N(); i++) count[i] = (Data::V(i) == Data::Fillval()) ? 0 : 1; + for(size_t i = 0; i < Data::N(); i++) count[i] = DataOk(this, i) ? 1 : 0; return; } if(Data::N() != d.N()) return; for(size_t i = 0; i < Data::N(); i++) - if(d(i) != Data::Fillval()) + if(DataOk(&d, i)) { - if(Data::V(i) == Data::Fillval()) + if(DataOk(this, i)) { - Data::V(i) = d(i); - count[i] = 1; + if constexpr(haveu) Data::U(i) += d.U(i); + if constexpr(havev) Data::V(i) += d.V(i); + count[i]++; } else { - Data::V(i) += d(i); - count[i]++; + if constexpr(haveu) Data::U(i) = d.U(i); + if constexpr(havev) Data::V(i) = d.V(i); + count[i] = 1; } } } @@ -86,7 +113,16 @@ template class Averager: public Data { if(!*this) return *this; for(size_t i = 0; i < Data::N(); i++) - if(count[i] != 0) Data::V(i) /= count[i]; + { + if constexpr(haveu) + { + if(count[i] != 0) Data::U(i) /= count[i]; + } + if constexpr(havev) + { + if(count[i] != 0) Data::V(i) /= count[i]; + } + } return *this; } }; diff --git a/include/layereddata.h b/include/layereddata.h index 0943721..99458ae 100644 --- a/include/layereddata.h +++ b/include/layereddata.h @@ -4,6 +4,7 @@ #include "gsw.h" #include "mdatetime.h" #include "simple2ddata.h" +#include "uvdata.h" #include #include #include @@ -148,6 +149,8 @@ class LayeredData Data Read(const MString& vname, const BaseParameters* ip, size_t i) const; + UVData ReadUV(const BaseParameters* ip, size_t i) const; + bool isOk() const { return nc.size() > 0; } explicit operator bool() const { return nc.size() > 0; } @@ -242,6 +245,8 @@ class LayeredData 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/odm.h b/include/odm.h index 21a95e1..f7c7a41 100644 --- a/include/odm.h +++ b/include/odm.h @@ -40,12 +40,9 @@ concept ReadSupported = requires(T t, const MString& vname, size_t i) { }; template -concept ActionTscSupported = requires(T t, michlib_internal::ParameterListEx& pars, const CLArgs& args, const BaseParameters* ip, size_t i) { +concept ReadUVPSupported = requires(T t, const BaseParameters* ip, size_t i) { { - t.template Parameters(pars, args).first - } -> std::convertible_to; - { - t.template Read(ip, i)(0) + t.template ReadUV(ip, i).U(0, 0) + t.template ReadUV(ip, i).V(0, 0) } -> std::convertible_to; }; @@ -163,4 +160,5 @@ class Data: public DataVariants MString ActionInfo(const CLArgs& args); MString ActionTsc(const CLArgs& args); + MString ActionUV(const CLArgs& args); }; diff --git a/include/uvdata.h b/include/uvdata.h new file mode 100644 index 0000000..4b16e39 --- /dev/null +++ b/include/uvdata.h @@ -0,0 +1,225 @@ +#pragma once +#include "geohelpers.h" +#include + +using michlib::M_PI; +using michlib::real; +using michlib::Sqrt; + +enum class Metric +{ + EUCLID, + SPHERIC +}; + +class UVData +{ + static constexpr real fillval = 1.0e10; + real x0 = 0.0, y0 = 0.0; + size_t nx = 0, ny = 0; + real xstep = 0.0, ystep = 0.0; + + std::vector u, v; + Metric metric; + + real D(real lon1, real lat1, real lon2, real lat2) const + { + switch(metric) + { + case(Metric::EUCLID): return michlib::Hypot(lon2 - lon1, lat2 - lat1); + // Spheric distance in km + case(Metric::SPHERIC): return michlib::GCD(lon1 * M_PI / 180.0, lat1 * M_PI / 180.0, lon2 * M_PI / 180.0, lat2 * M_PI / 180.0) * 6371.0; + } + return 0.0; + } + + // For derivatives + real Ud(size_t ix, size_t iy) const { return IsCoast(ix, iy) ? 0.0 : U(ix, iy); } + real Vd(size_t ix, size_t iy) const { return IsCoast(ix, iy) ? 0.0 : V(ix, iy); } + + real dUdX(size_t ix, size_t iy) const + { + if(IsCoast(ix, iy)) return fillval; + if(ix == 0) return (Ud(1, iy) - Ud(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy)); + if(ix == nx - 1) return (Ud(nx - 1, iy) - Ud(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy)); + return 0.5 * ((Ud(ix + 1, iy) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) + + (Ud(ix, iy) - Ud(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy))); + } + real dVdX(size_t ix, size_t iy) const + { + if(IsCoast(ix, iy)) return fillval; + if(ix == 0) return (Vd(1, iy) - Vd(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy)); + if(ix == nx - 1) return (Vd(nx - 1, iy) - Vd(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy)); + return 0.5 * ((Vd(ix + 1, iy) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) + + (Vd(ix, iy) - Vd(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy))); + } + real dUdY(size_t ix, size_t iy) const + { + if(IsCoast(ix, iy)) return fillval; + if(iy == 0) return (Ud(ix, 1) - Ud(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1)); + if(iy == ny - 1) return (Ud(ix, ny - 1) - Ud(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1)); + return 0.5 * ((Ud(ix, iy + 1) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) + + (Ud(ix, iy) - Ud(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy))); + } + real dVdY(size_t ix, size_t iy) const + { + if(IsCoast(ix, iy)) return fillval; + if(iy == 0) return (Vd(ix, 1) - Vd(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1)); + if(iy == ny - 1) return (Vd(ix, ny - 1) - Vd(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1)); + return 0.5 * ((Vd(ix, iy + 1) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) + + (Vd(ix, iy) - Vd(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy))); + } + + // For stationary points + real LonR(real x, [[maybe_unused]] real y) const { return x0 + x * xstep; } + real LatR([[maybe_unused]] real x, real y) const { return y0 + y * ystep; } + + public: + enum STPOINT + { + SADDLE = 0, + SACICFOCUS, + SKNOT, + UACICFOCUS, + UKNOT, + SCICFOCUS, + UCICFOCUS, + NOPOINT + }; + + struct StPoint + { + real x = 0.0, y = 0.0; + STPOINT type = NOPOINT; + }; + + UVData() = default; + UVData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_, Metric m_ = Metric::SPHERIC): + x0(x0_), y0(y0_), nx(nx_), ny(ny_), xstep(xs_), ystep(ys_), u(nx_ * ny_), v(nx_ * ny_), metric(m_) + { + } + + const real& U(size_t i) const { return u[i]; } + const real& V(size_t i) const { return v[i]; } + const real& U(size_t ix, size_t iy) const { return U(iy * nx + ix); } + const real& V(size_t ix, size_t iy) const { return V(iy * nx + ix); } + + real& U(size_t i) { return u[i]; } + real& V(size_t i) { return v[i]; } + real& U(size_t ix, size_t iy) { return U(iy * nx + ix); } + real& V(size_t ix, size_t iy) { return V(iy * nx + ix); } + + size_t N() const { return u.size(); } + size_t Nx() const { return nx; } + size_t Ny() const { return ny; } + + real Lon(size_t ix, [[maybe_unused]] size_t iy) const { return x0 + ix * xstep; } + 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); } + + explicit operator bool() const { return N() != 0; } + + bool IsCoast(size_t i) const { return U(i) == fillval || V(i) == fillval; } + bool IsCoast(size_t ix, size_t iy) const { return U(ix, iy) == fillval || V(ix, iy) == fillval; } + + static real Fillval() { return fillval; } + + real Div(size_t i) const { return Div(i % nx, i / nx); } + real Rot(size_t i) const { return Rot(i % nx, i / nx); } + real Div(size_t ix, size_t iy) const + { + if(!*this) return 0.0; + real ux = dUdX(ix, iy); + real vy = dVdY(ix, iy); + return (ux == fillval || vy == fillval) ? fillval : (ux + vy); + } + real Rot(size_t ix, size_t iy) const + { + if(!*this) return 0.0; + real vx = dVdX(ix, iy); + real uy = dUdY(ix, iy); + return (vx == fillval || uy == fillval) ? fillval : (vx - uy); + } + // Okubo-Weiss parameter + real OW(size_t i) const { return OW(i % nx, i / nx); } + real OW(size_t ix, size_t iy) const + { + if(!*this) return 0.0; + real ux = dUdX(ix, iy); + real uy = dUdY(ix, iy); + real vx = dVdX(ix, iy); + real vy = dVdY(ix, iy); + + if(ux == fillval || uy == fillval || vx == fillval || vy == fillval) return fillval; + + real sn = ux - vy; + real ss = vx + uy; + real w = vx - uy; + return sn * sn + ss * ss - w * w; + } + + auto StablePoints(size_t ix, size_t iy) const + { + std::vector points; + if(!*this) return points; + if(ix >= Nx() - 1 || iy >= Ny() - 1) return points; + if(IsCoast(ix, iy) || IsCoast(ix + 1, iy) || IsCoast(ix, iy + 1) || IsCoast(ix + 1, iy + 1)) return points; + + real au = U(ix, iy) - U(ix, iy + 1) - U(ix + 1, iy) + U(ix + 1, iy + 1); + real bu = U(ix, iy + 1) - U(ix, iy); + real cu = U(ix + 1, iy) - U(ix, iy); + real du = U(ix, iy); + + real av = V(ix, iy) - V(ix, iy + 1) - V(ix + 1, iy) + V(ix + 1, iy + 1); + real bv = V(ix, iy + 1) - V(ix, iy); + real cv = V(ix + 1, iy) - V(ix, iy); + real dv = V(ix, iy); + + real D = au * bv - av * bu; + real a = (au * av * cu - au * au * cv) / D; + real b = (au * bv * cu - au * bu * cv + au * av * du - au * au * dv) / D; + real c = (au * bv * du - au * bu * dv) / D; + real des = b * b - 4.0 * a * c; + if(des < 0) return points; + des = Sqrt(des); + + real p = av * cu - au * cv; + real q = av * bu - au * bv; + real r = av * du - au * dv; + real x1 = 0.5 * (-b + des) / a; + real x2 = 0.5 * (-b - des) / a; + real y1 = -(p * x1 + r) / q; + real y2 = -(p * x2 + r) / q; + + auto PointType = [au = au, av = av, bu = bu, bv = bv, cu = cu, cv = cv](real x, real y) -> auto + { + real a = 1.0; + real b = -y * au - x * av - bv - cu; + real c = y * (au * bv - av * bu) + x * (av * cu - au * cv) + bv * cu - bu * cv; + real des = b * b - 4.0 * a * c; + if(des > 0.0) + { + if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0 && c > 0.0) + return UKNOT; + else if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0) + return SADDLE; + else + return SKNOT; + } + else + { + bool acyclcrit = (av * y + cv > 0.0); + if(b < 0.0) + return acyclcrit ? SACICFOCUS : SCICFOCUS; + else + return acyclcrit ? UACICFOCUS : UCICFOCUS; + } + }; + + if(x1 >= 0.0 && x1 < 1.0 && y1 >= 0.0 && y1 < 1.0) points.emplace_back(LonR(x1 + ix, y1 + iy), LatR(x1 + ix, y1 + iy), PointType(x1, y1)); + if(x2 >= 0.0 && x2 < 1.0 && y2 >= 0.0 && y2 < 1.0) points.emplace_back(LonR(x2 + ix, y2 + iy), LatR(x2 + ix, y2 + iy), PointType(x2, y2)); + return points; + } +}; diff --git a/src/actiontsc.cpp b/src/actiontsc.cpp index afad0b6..14a7abc 100644 --- a/src/actiontsc.cpp +++ b/src/actiontsc.cpp @@ -96,6 +96,9 @@ MString Data::ActionTsc(const CLArgs& args) if(!name.Exist()) name = "out.bin"; fw.Create(name, 3); + fw.SetColumnName(1, "Longitude"); + fw.SetColumnName(2, "Latitude"); + fw.SetColumnName(3, vname); fw.SetParameters(pars); for(size_t i = 0; i < data.N(); i++) { diff --git a/src/actionuv.cpp b/src/actionuv.cpp new file mode 100644 index 0000000..8c3fb0f --- /dev/null +++ b/src/actionuv.cpp @@ -0,0 +1,169 @@ +#define MICHLIB_NOSOURCE +#include "GPL.h" +#include "odm.h" + +MString Data::ActionUV(const CLArgs& args) +{ + if(args.contains("time") && (args.contains("timeb") || args.contains("timee"))) return "Time must be set via time parameter or timeb and timee parameter but not via both"; + if(!(args.contains("time") || (args.contains("timeb") && args.contains("timee")))) return "Time must be set via time parameter or timeb and timee parameter"; + + michlib_internal::ParameterListEx pars; + pars.UsePrefix(""); + + pars.SetParameter("source", args.at("source")); + TIndex tindexes; + + if(args.contains("time")) + { + tindexes = GetTIndexes(args.at("time")); // Regular expression or single date + if(tindexes.size() == 0) return "There are no times matching the regular expression"; + if(tindexes.size() == 1) + pars.SetParameter("time", Time(tindexes[0]).ToString()); + else + pars.SetParameter("timeregex", args.at("time")); + } + else + { + MDateTime tb(args.at("timeb")), te(args.at("timee")); + tindexes = GetTIndexes(tb, te); + if(tindexes.size() == 0) return "There are no times between " + tb.ToString() + " and " + te.ToString(); + pars.SetParameter("timeb", tb.ToString()); + pars.SetParameter("timee", te.ToString()); + } + + std::unique_ptr sourcepars; + { + auto ppar = std::visit( + [&pars = pars, &args = args](const auto& arg) -> std::pair + { + using T = std::decay_t; + if constexpr(ParametersSupported) + return arg.Parameters(pars, args); + else + return {nullptr, ""}; + }, + *this); + if(ppar.second.Exist()) return ppar.second; + sourcepars.reset(ppar.first); + } + + auto p = sourcepars.get(); + size_t ind; + auto read = [p = p, &ind = ind](const auto& arg) -> UVData + { + using T = std::decay_t; + if constexpr(ReadUVPSupported) + return arg.ReadUV(p, ind); + else + return UVData(); + }; + + UVData data; + if(tindexes.size() == 1) + { + ind = tindexes[0]; + michlib::message("Time: " + Time(ind).ToTString()); + data = std::visit(read, *this); + } + else + { + Averager out; + bool ok = true; + for(size_t i = 0; i < tindexes.size(); i++) + { + if(!ok) break; + ind = tindexes[i]; + michlib::message("Time: " + Time(ind).ToTString()); + auto dat = std::visit(read, *this); + if(dat) + out.Add(dat); + else + ok = false; + } + if(ok) data = out.Div(); + } + if(!data) return "Can't read data"; + + // Main file + MString name = args.contains("out") ? args.at("out") : ""; + if(name.Exist()) + { + BFileW fw; + fw.Create(name, 7); + fw.SetColumnName(1, "Longitude"); + fw.SetColumnName(2, "Latitude"); + fw.SetColumnName(3, "u, cm/s"); + fw.SetColumnName(4, "v, cm/s"); + fw.SetColumnName(5, "Div, (cm/s)/km"); + fw.SetColumnName(6, "Rot, (cm/s)/km"); + fw.SetColumnName(7, "Okubo-Weiss parameter, (cm^2/s^2)/km^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.Finalize(); + fw.Close(); + } + + // Filtered vectors file + name = args.contains("velout") ? args.at("velout") : ""; + if(name.Exist()) + { + 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; + + BFileW vel; + vel.Create(name, 4); + vel.SetColumnName(1, "Longitude"); + vel.SetColumnName(2, "Latitude"); + vel.SetColumnName(3, "u, cm/s"); + vel.SetColumnName(4, "v, cm/s"); + + for(size_t ix = shiftx; ix < data.Nx(); ix += skipx) + for(size_t iy = shifty; iy < data.Ny(); iy += skipy) + { + vel.Write(data.Lon(ix, iy)); + vel.Write(data.Lat(ix, iy)); + vel.Write(data.U(ix, iy) == data.Fillval() ? NAN : data.U(ix, iy)); + vel.Write(data.V(ix, iy) == data.Fillval() ? NAN : data.V(ix, iy)); + } + vel.Finalize(); + vel.Close(); + } + + // Stationary points + name = args.contains("stpout") ? args.at("stpout") : ""; + if(name.Exist()) + { + BFileW stp; + stp.Create(name, 3); + stp.SetColumnName(1, "Longitude"); + stp.SetColumnName(2, "Latitude"); + stp.SetColumnName(3, "Stability (0 - saddle, 1 - st. anticicl. focus, 2 - st. knot, 3 - unst. anticicl. focus, 4 - unst. knot, 5 - st. cicl. focus, 6 - unst. cicl. focus)"); + + for(size_t ix = 0; ix < data.Nx() - 1; ix++) + for(size_t iy = 0; iy < data.Ny() - 1; iy++) + { + auto points = data.StablePoints(ix, iy); + for(size_t i = 0; i < points.size(); i++) + { + stp.Write(points[i].x); + stp.Write(points[i].y); + stp.Write(michlib::int_cast(points[i].type)); + } + } + stp.Finalize(); + stp.Close(); + } + + return ""; +} diff --git a/src/layereddata.cpp b/src/layereddata.cpp index 559a356..02d0e91 100644 --- a/src/layereddata.cpp +++ b/src/layereddata.cpp @@ -280,6 +280,28 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters* return Data(); } +UVData LayeredData::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), lonstep, latstep}; + 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; +} + template LayeredData::Data LayeredData::ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct LayeredData::Parameters* p) const { real unitmul = 1.0; diff --git a/src/odm.cpp b/src/odm.cpp index 1dd7105..62dbfa7 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. Default: tsc."); + 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"); message(" Keys for source=NEMO"); @@ -12,7 +12,7 @@ inline void Usage(const MString& arg0) message(" Keys for source=HYCOM"); message(" dataset. Can be Forecast, Hindcast or Reanalysis. Default: Forecast"); message(" Keys for action=tsc. Get temperature, salinity, chlorofill from dataset."); - message(" source. Required. May be: NEMO"); + message(" source. Required. May be: NEMO, HYCOM"); message(" var. Required. May be: 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"); @@ -25,6 +25,17 @@ inline void Usage(const MString& arg0) message(" dataset. Can be Forecast, Hindcast or Reanalysis. Default: Forecast"); message(" layer and/or depth. Layer or depth of HYCOM dataset. If depth is specified, layer is ignored. Both ignored, if var=mld or var=ssh. Default: layer=0"); message(" lonb, lone, latb, late. Required. Region of interest"); + message(" Keys for action=uv. Get velocity field and its derivatives."); + message(" source. Required. May be: NEMO, HYCOM"); + 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."); + message(" velout. Output file for sparse velocity field. If present, parameters shiftx, shifty, skipx, skipy controls the sparsing."); + message(" shiftx, shifty. Shift of the sparsed grid. Used only if velout parameter is present. Default: 0"); + 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 source=NEMO. See section action=tsc."); + message(" Keys for source=HYCOM. See section action=tsc."); } int main(int argc, char** argv) @@ -53,6 +64,11 @@ int main(int argc, char** argv) ret = data.ActionInfo(args); else if(action == "tsc") ret = data.ActionTsc(args); + else if(action == "uv") + { + args["var"] = "U"; + ret = data.ActionUV(args); + } else { errmessage("Unknown action " + action);