You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

335 lines
12 KiB

#pragma once
#include "BFileW.h"
#include "actiondep.h"
#include "ncfilew.h"
#include <utility>
2 years ago
using michlib::BFileW;
class UVMethods
{
protected:
template<class D> class SparserHolder
{
protected:
const D& data;
size_t shiftx, shifty, skipx, skipy;
SparserHolder(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): data(d), shiftx(shx), shifty(shy), skipx(sx), skipy(sy) {}
size_t I2Ix(size_t i) const { return i * skipx + shiftx; }
size_t I2Iy(size_t i) const { return i * skipy + shifty; }
public:
size_t Nx() const
{
size_t n = data.Nx() - shiftx;
return n / skipx + ((n % skipx > 0) ? 1 : 0);
}
size_t Ny() const
{
size_t n = data.Ny() - shifty;
return n / skipy + ((n % skipy > 0) ? 1 : 0);
}
static constexpr auto Fillval() { return D::Fillval(); }
real Lon(size_t ix, size_t iy) const { return data.Lon(I2Ix(ix), I2Iy(iy)); }
real Lat(size_t ix, size_t iy) const { return data.Lat(I2Ix(ix), I2Iy(iy)); }
real U(size_t ix, size_t iy) const { return data.U(I2Ix(ix), I2Iy(iy)); }
real V(size_t ix, size_t iy) const { return data.V(I2Ix(ix), I2Iy(iy)); }
bool IsFill(size_t ix, size_t iy) const { return data.IsFill(I2Ix(ix), I2Iy(iy)); }
};
template<class D, bool b> class Sparser2DGeoRectArray: public SparserHolder<D>
{
protected:
Sparser2DGeoRectArray(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): SparserHolder<D>(d, shx, shy, sx, sy) {}
};
template<class D> class Sparser2DGeoRectArray<D, true>: public SparserHolder<D>
{
protected:
using SparserHolder<D>::data;
using SparserHolder<D>::I2Ix;
using SparserHolder<D>::I2Iy;
Sparser2DGeoRectArray(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): SparserHolder<D>(d, shx, shy, sx, sy) {}
public:
real Ix2Lon(size_t ix) const { return data.Ix2Lon(I2Ix(ix)); }
real Iy2Lat(size_t iy) const { return data.Iy2Lat(I2Iy(iy)); }
};
template<class D> class Sparser: public Sparser2DGeoRectArray<D, ReadIs2DGeoRectArray<D>>
{
public:
Sparser(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): Sparser2DGeoRectArray<D, ReadIs2DGeoRectArray<D>>(d, shx, shy, sx, sy) {}
};
class StPoints
{
std::vector<real> x, y;
std::vector<michlib::uint1> t;
bool lonlat;
NCFileWBase nc;
bool tdep;
size_t N() const { return t.size(); }
public:
StPoints(bool ll): lonlat(ll) {}
void Add(const std::vector<UVDataStorage::StPoint>& points)
{
for(const auto& p: points)
{
x.push_back(p.x);
y.push_back(p.y);
t.push_back(michlib::int_cast<michlib::uint1>(p.type));
}
}
void Reset()
{
x.clear();
y.clear();
t.clear();
}
void WriteBinBile(const MString& name, const michlib_internal::ParameterListEx& pars) 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(); }
};
};
ADD_ACTION(UV, uv, ReadPSupported<Source> || ReadSupported<Source>, UVMethods);
2 years ago
template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
{
auto [reg, regerr] = GetRegion<D>(args);
if(regerr.Exist()) return regerr;
auto resop = ds.Open(args);
if(resop.Exist()) return "Can't open source: " + resop;
2 years ago
michlib_internal::ParameterListEx pars;
pars.UsePrefix("");
pars.SetParameter("source", args.at("source"));
pars.SetParameter("history", args.at("_cmdline"));
2 years ago
bool mode = args.contains("geostrophic");
auto [tindexes, err] = GetTIndexes(ds, args, pars);
if(err.Exist()) return err;
2 years ago
bool average = args.contains("average");
int compress = 3;
if(args.contains("compress")) compress = args.at("compress").ToInt();
2 years ago
std::unique_ptr<const BaseParameters> sourcepars;
if constexpr(ParametersSupported<D>)
2 years ago
{
if constexpr(ParametersRequiredRegion<D>)
{
auto [p, err] = ds.Parameters(pars, args, reg);
if(err.Exist()) return err;
sourcepars.reset(p);
}
else
{
auto [p, err] = ds.Parameters(pars, args);
if(err.Exist()) return err;
sourcepars.reset(p);
}
2 years ago
}
auto p = sourcepars.get();
2 years ago
TimeData tdata(ds, tindexes);
2 years ago
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(namevel) == "nc" ? "nc" : "bin");
MString namestp = args.contains("stpout") ? args.at("stpout") : "";
MString outfmtstp = args.contains("stpoutformat") ? args.at("stpoutformat") : (GetExt(namestp) == "nc" ? "nc" : "bin");
size_t shiftx = args.contains("shiftx") ? args.at("shiftx").ToInteger<size_t>() : 0;
size_t shifty = args.contains("shifty") ? args.at("shifty").ToInteger<size_t>() : 0;
size_t skipx = args.contains("skipx") ? args.at("skipx").ToInteger<size_t>() : 1;
size_t skipy = args.contains("skipy") ? args.at("skipy").ToInteger<size_t>() : 1;
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";
bool headwrited = false;
MString velunit, distunit;
NCFileW fw, fwfilt;
StPoints stp(true);
2 years ago
for(size_t it = 0; it < tindexes.size(); it++)
2 years ago
{
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";
}
2 years ago
// Main file
if(name.Exist())
{
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++)
{
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")
{
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, <u^2+v^2>-<u>^2-<v>^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;
}
// Filtered file
if(namevel.Exist())
{
Sparser<decltype(data)> sdata(data, shiftx, shifty, skipx, skipy);
if(outfmtvel == "bin")
{
BFileW vel;
vel.Create(namevel, 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, namevel, 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;
}
// Stationary points
if(namestp.Exist())
{
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(namestp, pars);
else if(outfmtstp == "nc" || outfmtstp == "netcdf")
{
MString err;
if(!err.Exist() && !headwrited) err = stp.CreateNcFile(namestp, pars, args.at("_cmdline"), compress, tdata, !average);
if(!err.Exist()) err = stp.WriteNcFile(it);
if(err.Exist()) return err;
}
else
return "Unknown format: " + outfmt;
}
2 years ago
if(average) break;
headwrited = true;
}
fw.Close();
fwfilt.Close();
stp.CloseNcFile();
2 years ago
return "";
};