Compare commits

...

13 Commits

Author SHA1 Message Date
Michael Uleysky 18b1f939b3 New class for access to vylet data 9 months ago
Michael Uleysky 2e0f396497 A new class Rect2DData for data along meridians and parallels, but not necessarily with uniform spacing. 9 months ago
Michael Uleysky 9f4d7094dd Sources may specify default variable lists for the tsc module. 10 months ago
Michael Uleysky df85e7ca51 Reading and caching multiple variables from the source. 10 months ago
Michael Uleysky c85dc28858 Removed copy constructor from BaseData to prevent unnecessary copying of data 10 months ago
Michael Uleysky 3b40409fb2 Sources can now set the standart_name, long_name, and comment attributes on read data. If the standart_name and long_name attributes are not set, the tsc module will attempt to set them according to the variable names. 10 months ago
Michael Uleysky ed1282413f Added support for netcdf output for the uv action 10 months ago
Michael Uleysky 1d9268b32d Using the new NCFileW class to write netcdf files with the tsc module 10 months ago
Michael Uleysky 3b1d5c3157 New NCFileW class to simplify writing netcdf files 10 months ago
Michael Uleysky fd7d7f8a5f Interfaces of UVData classes made more consistent with interfaces of BaseData classes 10 months ago
Michael Uleysky 8a1cf03804 Improved ReadType-based concepts 10 months ago
Michael Uleysky e25f843e7e Autodetect format by output file extension in the tsc action 10 months ago
Michael Uleysky 20824cdaab Added unit attribute in the netcdf output file. 10 months ago
  1. 191
      actions/actiontsc.h
  2. 88
      actions/actionuv.cpp
  3. 277
      actions/actionuv.h
  4. 1
      include/ParseArgs.h
  5. 85
      include/actiondep.h
  6. 18
      include/basedata.h
  7. 2
      include/layereddata.h
  8. 201
      include/ncfilew.h
  9. 35
      include/simple2ddata.h
  10. 61
      include/traits.h
  11. 22
      include/uvdata.h
  12. 43
      sources/AVISOLOCAL.cpp
  13. 2
      sources/AVISOLOCAL.h
  14. 32
      sources/BINFILE.cpp
  15. 2
      sources/BINFILE.h
  16. 12
      sources/MODISBINLOCAL.cpp
  17. 11
      sources/MODISBINLOCAL.h
  18. 526
      sources/VYLET.cpp
  19. 145
      sources/VYLET.h
  20. 7
      src/ParseArgs.cpp
  21. 105
      src/layereddata.cpp
  22. 105
      src/ncfilew.cpp

191
actions/actiontsc.h

@ -1,6 +1,7 @@
#pragma once
#include "BFileW.h"
#include "actiondep.h"
#include "ncfilew.h"
#include "ncfuncs.h"
#include <memory>
@ -23,10 +24,35 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
auto [tindexes, err] = GetTIndexes(ds, args, pars);
if(err.Exist()) return err;
if(!args.contains("var")) return "Variable not specified";
MString vname = args.at("var");
if(!ds.CheckVar(vname)) return "Variable " + vname + " not exists in this dataset";
pars.SetParameter("variable", vname);
MString varstring;
if(args.contains("var"))
varstring = args.at("var");
else
{
if(args.contains("vars"))
varstring = args.at("vars");
else
{
if constexpr(HasDefVars<D>)
varstring = ds.DefaultVars();
else
return "Variables not specified";
}
}
auto vlist = michlib::Split_on_words(varstring, " ,", false);
std::vector<MString> vnames{std::make_move_iterator(std::begin(vlist)), std::make_move_iterator(std::end(vlist))};
{
std::set<MString> used;
for(const auto& vname: vnames)
{
if(used.contains(vname)) return "Duplicate variable " + vname + " in list " + varstring;
if(!ds.CheckVar(vname)) return "Variable " + vname + " not exists in this dataset";
used.insert(vname);
}
}
pars.SetParameter("variables", varstring);
std::unique_ptr<const BaseParameters> sourcepars;
if constexpr(ParametersSupported<D>)
@ -46,27 +72,34 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
}
auto p = sourcepars.get();
auto data = Read(ds, vname, p, tindexes);
if(!data) return "Can't read data";
if(!data.Unit().Exist()) michlib::errmessage("Unknown measurement unit!");
auto data = Read(ds, vnames, p, tindexes);
if(data.size() != vnames.size()) return "Can't read data";
MString outfmt = args.contains("format") ? args.at("format") : "bin";
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].SetStandartName(NCFuncs::Name2LongName(vnames[i]));
}
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")
{
BFileW fw;
MString name = args.contains("out") ? args.at("out") : "out.bin";
BFileW fw;
fw.Create(name, 3);
fw.SetColumnName(1, "Longitude");
fw.SetColumnName(2, "Latitude");
fw.SetColumnName(3, vname + ", " + (data.Unit().Exist() ? data.Unit() : "unknown"));
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 i = 0; i < data.N(); i++)
for(size_t r = 0; r < data[0].N(); r++)
{
fw.Write(data.Lon(i));
fw.Write(data.Lat(i));
fw.Write(data.IsFill(i) ? NAN : data(i));
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();
@ -75,122 +108,20 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
if(outfmt == "nc" || outfmt == "netcdf")
{
MString name = args.contains("out") ? args.at("out") : "out.nc";
int ret;
int ncid, dimid, dimidxy[2];
int lonid, latid, vid;
MString text;
auto fill = std::numeric_limits<float>::max();
float fval;
int compress = 3;
NCFileW fw;
MString err;
if(args.contains("compress")) compress = args.at("compress").ToInt();
ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &ncid);
if(ret != NC_NOERR) return "Can't create netcdf file: " + name;
ret = nc_put_att_text(ncid, NC_GLOBAL, "history", args.at("_cmdline").Len() + 1, args.at("_cmdline").Buf());
if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
if constexpr(ReadIs2DGeoRectArray<D>)
{
ret = nc_def_dim(ncid, "longitude", data.Nx(), &dimidxy[0]);
if(ret != NC_NOERR) return "Can't create x-dimension in the netcdf file";
ret = nc_def_dim(ncid, "latitude", data.Ny(), &dimidxy[1]);
if(ret != NC_NOERR) return "Can't create y-dimension in the netcdf file";
ret = nc_def_var(ncid, "longitude", NC_FLOAT, 1, &dimidxy[0], &lonid);
if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file";
ret = nc_def_var(ncid, "latitude", NC_FLOAT, 1, &dimidxy[1], &latid);
if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file";
ret = nc_def_var(ncid, vname.Buf(), NC_FLOAT, 2, dimidxy, &vid);
if(ret != NC_NOERR) return "Can't create " + vname + " variable in the netcdf file";
}
else
{
ret = nc_def_dim(ncid, "i", data.N(), &dimid);
if(ret != NC_NOERR) return "Can't create dimension in the netcdf file";
ret = nc_def_var(ncid, "longitude", NC_FLOAT, 1, &dimid, &lonid);
if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file";
ret = nc_def_var(ncid, "latitude", NC_FLOAT, 1, &dimid, &latid);
if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file";
ret = nc_def_var(ncid, vname.Buf(), NC_FLOAT, 1, &dimid, &vid);
if(ret != NC_NOERR) return "Can't create " + vname + " variable in the netcdf file";
}
ret = nc_def_var_deflate(ncid, lonid, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for longitude variable in the netcdf file";
ret = nc_def_var_deflate(ncid, latid, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for latitude variable in the netcdf file";
ret = nc_def_var_deflate(ncid, vid, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for " + vname + " variable in the netcdf file";
text = "longitude";
ret = nc_put_att_text(ncid, lonid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of longitude variable in the netcdf file";
text = "latitude";
ret = nc_put_att_text(ncid, latid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of latitude variable in the netcdf file";
text = NCFuncs::Name2StName(vname);
if(text.Exist())
{
ret = nc_put_att_text(ncid, vid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of " + vname + " variable in the netcdf file";
}
else
{
text = NCFuncs::Name2LongName(vname);
ret = nc_put_att_text(ncid, vid, "long_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of " + vname + " variable in the netcdf file";
}
ret = nc_put_att_float(ncid, vid, "_FillValue", NC_FLOAT, 1, &fill);
if(ret != NC_NOERR) return "Can't write _FillValue attribute of " + vname + " variable in the netcdf file";
ret = nc_enddef(ncid);
if(ret != NC_NOERR) return "Can't finish definition of the netcdf file";
if constexpr(ReadIs2DGeoRectArray<D>)
{
size_t i[2];
for(i[0] = 0; i[0] < data.Nx(); i[0]++)
{
fval = data.Ix2Lon(i[0]);
ret = nc_put_var1_float(ncid, lonid, &i[0], &fval);
if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file: " + MString(i[0]);
}
for(i[1] = 0; i[1] < data.Ny(); i[1]++)
{
fval = data.Iy2Lat(i[1]);
ret = nc_put_var1_float(ncid, latid, &i[1], &fval);
if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file: " + MString(i[1]);
}
for(i[0] = 0; i[0] < data.Nx(); i[0]++)
for(i[1] = 0; i[1] < data.Ny(); i[1]++)
{
fval = data.IsFill(i[0], i[1]) ? fill : data(i[0], i[1]);
ret = nc_put_var1_float(ncid, vid, i, &fval);
if(ret != NC_NOERR) return "Can't write " + vname + " variable in the netcdf file: " + MString(i[0]) + ", " + i[1];
}
}
else
{
for(size_t i = 0; i < data.N(); i++)
{
fval = data.Lon(i);
ret = nc_put_var1_float(ncid, lonid, &i, &fval);
if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file: " + MString(i);
fval = data.Lat(i);
ret = nc_put_var1_float(ncid, latid, &i, &fval);
if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file: " + MString(i);
fval = data.IsFill(i) ? fill : data(i);
ret = nc_put_var1_float(ncid, vid, &i, &fval);
if(ret != NC_NOERR) return "Can't write " + vname + " variable in the netcdf file: " + MString(i);
}
}
ret = nc_close(ncid);
if(ret != NC_NOERR) return "Can't close netcdf file";
if(!err.Exist()) err = fw.Create(data[0], name, args.at("_cmdline"), compress);
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 "";
}

88
actions/actionuv.cpp

@ -0,0 +1,88 @@
#define MICHLIB_NOSOURCE
#include "actionuv.h"
void UVMethods::StPoints::WriteBinBile(const MString& name) const
{
BFileW stp;
stp.Create(name, 3);
stp.SetColumnName(1, lonlat ? "Longitude" : "x");
stp.SetColumnName(2, lonlat ? "Latitude" : "x");
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 i = 0; i < N(); i++)
{
stp.Write(x[i]);
stp.Write(y[i]);
stp.Write(michlib::int_cast<size_t>(t[i]));
}
stp.Finalize();
stp.Close();
}
MString UVMethods::StPoints::WriteNcFile(const MString& name, const MString& history, int comp) const
{
int ret;
int ncid;
int dimid;
int xid, yid, tid;
MString text;
ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &ncid);
if(ret != NC_NOERR) return "Can't create netcdf file: " + name;
ret = nc_put_att_text(ncid, NC_GLOBAL, "history", history.Len() + 1, history.Buf());
if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
ret = nc_def_dim(ncid, "i", N(), &dimid);
if(ret != NC_NOERR) return "Can't create dimension in the netcdf file";
ret = nc_def_var(ncid, lonlat ? "longitude" : "x", NC_FLOAT, 1, &dimid, &xid);
if(ret != NC_NOERR) return "Can't create x variable in the netcdf file";
ret = nc_def_var(ncid, lonlat ? "latitude" : "y", NC_FLOAT, 1, &dimid, &yid);
if(ret != NC_NOERR) return "Can't create y variable in the netcdf file";
if(lonlat)
{
text = "longitude";
ret = nc_put_att_text(ncid, xid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of longitude variable in the netcdf file";
text = "latitude";
ret = nc_put_att_text(ncid, yid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of latitude variable in the netcdf file";
}
ret = nc_def_var_deflate(ncid, xid, 1, 1, comp);
if(ret != NC_NOERR) return "Can't set deflate parameters for x variable in the netcdf file";
ret = nc_def_var_deflate(ncid, yid, 1, 1, comp);
if(ret != NC_NOERR) return "Can't set deflate parameters for y variable in the netcdf file";
ret = nc_def_var(ncid, "type", NC_UBYTE, 1, &dimid, &tid);
if(ret != NC_NOERR) return "Can't create type variable in the netcdf file";
ret = nc_def_var_deflate(ncid, tid, 1, 1, comp);
if(ret != NC_NOERR) return "Can't set deflate parameters for type variable in the netcdf file";
text = "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";
ret = nc_put_att_text(ncid, tid, "long_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of type variable in the netcdf file";
ret = nc_enddef(ncid);
if(ret != NC_NOERR) return "Can't finish definition of the netcdf file";
const size_t i = 0, c = N();
float buf[c];
for(size_t ix = 0; ix < c; ix++) buf[ix] = x[ix];
ret = nc_put_vara_float(ncid, xid, &i, &c, buf);
if(ret != NC_NOERR) return "Can't write x variable in the netcdf file";
for(size_t ix = 0; ix < c; ix++) buf[ix] = y[ix];
ret = nc_put_vara_float(ncid, yid, &i, &c, buf);
if(ret != NC_NOERR) return "Can't write y variable in the netcdf file";
ret = nc_put_vara_ubyte(ncid, tid, &i, &c, t.data());
if(ret != NC_NOERR) return "Can't write type variable in the netcdf file";
ret = nc_close(ncid);
if(ret != NC_NOERR) return "Can't close netcdf file";
return "";
}

277
actions/actionuv.h

@ -1,11 +1,98 @@
#pragma once
#include "BFileW.h"
#include "actiondep.h"
#include <memory>
#include "ncfilew.h"
#include <utility>
using michlib::BFileW;
ADD_ACTION(UV, uv, ReadPSupported<Source> || ReadSupported<Source>);
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;
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 WriteBinBile(const MString& name) const;
MString WriteNcFile(const MString& name, const MString& history, int comp) const;
};
};
ADD_ACTION(UV, uv, ReadPSupported<Source> || ReadSupported<Source>, UVMethods);
template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
{
@ -44,39 +131,81 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
if(!data) return "Can't read data";
// Main file
MString name = args.contains("out") ? args.at("out") : "";
MString name = args.contains("out") ? args.at("out") : "";
MString outfmt = args.contains("outformat") ? args.at("outformat") : (GetExt(name) == "nc" ? "nc" : "bin");
MString u = data.Unit().Exist() ? data.Unit() : "unknown", d = data.DUnit().Exist() ? data.DUnit() : "unknown";
if(name.Exist())
{
BFileW fw;
fw.Create(name, 9);
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.SetColumnName(8, "Kinetic energy, cm^2/s^2");
fw.SetColumnName(9, "Eddy kinetic energy, cm^2/s^2");
fw.SetParameters(pars);
for(size_t i = 0; i < data.N(); i++)
if(outfmt == "bin")
{
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))));
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();
}
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, args.at("_cmdline"), compress);
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, <u^2+v^2>-<u>^2-<v>^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;
fw.Close();
}
else
return "Unknown format: " + outfmt;
}
// Filtered vectors file
name = args.contains("velout") ? args.at("velout") : "";
name = args.contains("velout") ? args.at("velout") : "";
outfmt = args.contains("veloutformat") ? args.at("veloutformat") : (GetExt(name) == "nc" ? "nc" : "bin");
if(name.Exist())
{
size_t shiftx = args.contains("shiftx") ? args.at("shiftx").ToInteger<size_t>() : 0;
@ -84,48 +213,70 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
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;
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");
Sparser<decltype(data)> sdata(data, shiftx, shifty, skipx, skipy);
if(outfmt == "bin")
{
BFileW vel;
vel.Create(name, 4);
vel.SetColumnName(1, "Longitude");
vel.SetColumnName(2, "Latitude");
vel.SetColumnName(3, "u, " + u);
vel.SetColumnName(4, "v, " + u);
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();
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(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, args.at("_cmdline"), compress);
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();
}
else
return "Unknown format: " + outfmt;
}
// Stationary points
name = args.contains("stpout") ? args.at("stpout") : "";
name = args.contains("stpout") ? args.at("stpout") : "";
outfmt = args.contains("stpoutformat") ? args.at("stpoutformat") : (GetExt(name) == "nc" ? "nc" : "bin");
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)");
StPoints p(true);
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<size_t>(points[i].type));
}
}
stp.Finalize();
stp.Close();
for(size_t iy = 0; iy < data.Ny() - 1; iy++) p.Add(data.StablePoints(ix, iy));
if(outfmt == "bin")
p.WriteBinBile(name);
else if(outfmt == "nc" || outfmt == "netcdf")
{
int compress = args.contains("compress") ? args.at("compress").ToInt() : 3;
MString err = p.WriteNcFile(name, args.at("_cmdline"), compress);
if(err.Exist()) return err;
}
else
return "Unknown format: " + outfmt;
}
return "";

1
include/ParseArgs.h

@ -6,3 +6,4 @@ using michlib::SList;
using CLArgs = std::map<MString, MString>;
CLArgs ParseArgs(int argc, char** argv);
MString GetExt(const MString& fname);

85
include/actiondep.h

@ -105,41 +105,50 @@ template<class D> std::pair<TIndex, MString> GetTIndexes(const D& data, const CL
return {tindexes, ""};
}
template<class D> ReadType<D> Read(const D& data, const MString& vname, const BaseParameters* p, const TIndex& tindex)
template<class D> std::vector<ReadType<D>> Read(const D& data, const std::vector<MString>& vnames, const BaseParameters* p, const TIndex& tindex)
{
using RT = ReadType<D>;
size_t ind;
size_t ind;
std::vector<RT> out;
if(tindex.size() == 1)
{
ind = tindex[0];
michlib::message("Time: " + data.Time(ind).ToTString());
if constexpr(ReadPSupported<D>)
return data.Read(vname, p, ind);
else if constexpr(ReadSupported<D>)
return data.Read(vname, ind);
std::map<MString, RT> cache;
for(const auto& vname: vnames)
{
bool res;
if constexpr(ReadPSupported<D>)
res = data.Read(vname, cache, p, ind);
else if constexpr(ReadSupported<D>)
res = data.Read(vname, cache, ind);
if(!res) return out;
}
for(size_t i = 0; i < vnames.size(); i++) out.emplace_back(std::move(cache[vnames[i]]));
}
else
{
Averager<RT> out;
bool ok = true;
std::vector<Averager<RT>> aver(vnames.size());
for(size_t i = 0; i < tindex.size(); i++)
{
if(!ok) break;
ind = tindex[i];
michlib::message("Time: " + data.Time(ind).ToTString());
RT dat;
if constexpr(ReadPSupported<D>)
dat = data.Read(vname, p, ind);
else if constexpr(ReadSupported<D>)
dat = data.Read(vname, ind);
if(dat)
out.Add(dat);
else
ok = false;
std::map<MString, RT> cache;
for(const auto& vname: vnames)
{
bool res;
if constexpr(ReadPSupported<D>)
res = data.Read(vname, cache, p, ind);
else if constexpr(ReadSupported<D>)
res = data.Read(vname, cache, ind);
if(!res) return out;
}
for(size_t i = 0; i < vnames.size(); i++) aver[i].Add(std::move(cache[vnames[i]]));
}
if(ok) return out.Div();
for(size_t i = 0; i < vnames.size(); i++) out.emplace_back(std::move(aver[i].Div()));
return out;
}
return RT();
return out;
}
template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, size_t ind)
@ -147,18 +156,15 @@ template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters
using RT = ReadType<D>;
using UV = UVData<RT>;
michlib::message("Time: " + data.Time(ind).ToTString());
RT u, v;
std::map<MString, RT> cache;
bool res = false;
if constexpr(ReadPSupported<D>)
{
u = data.Read("u", p, ind);
v = data.Read("v", p, ind);
}
res = data.Read("u", cache, p, ind) && data.Read("v", cache, p, ind);
else if constexpr(ReadSupported<D>)
{
u = data.Read("u", ind);
v = data.Read("v", ind);
}
return UV(u, v);
res = data.Read("u", cache, ind) && data.Read("v", cache, ind);
if(!res) return UV();
return UV(cache.at("u"), cache.at("v"));
}
template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex)
@ -178,20 +184,17 @@ template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters
if(!ok) break;
ind = tindex[i];
michlib::message("Time: " + data.Time(ind).ToTString());
RT u, v;
std::map<MString, RT> cache;
bool res = false;
if constexpr(ReadPSupported<D>)
{
u = data.Read("u", p, ind);
v = data.Read("v", p, ind);
}
res = data.Read("u", cache, p, ind) && data.Read("v", cache, p, ind);
else if constexpr(ReadSupported<D>)
{
u = data.Read("u", ind);
v = data.Read("v", ind);
}
UV dat(u, v);
res = data.Read("u", cache, ind) && data.Read("v", cache, ind);
if(!res) return UV();
UV dat(cache.at("u"), cache.at("v"));
if(dat)
out.Add(dat);
out.Add(std::move(dat));
else
ok = false;
}

18
include/basedata.h

@ -20,12 +20,16 @@ class BaseData
protected:
static constexpr real fillval = 1.0e10;
std::vector<real> data;
MString unit;
MString unit, stname, lname, comment;
BaseData(size_t n, MString&& unit_ = ""): data(n), unit(std::move(unit_)) {}
public:
BaseData() = default;
BaseData() = default;
BaseData(const BaseData&) = delete;
BaseData(BaseData&&) = default;
BaseData& operator=(BaseData&&) = default;
const real& V(size_t i) const { return data[i]; }
real& V(size_t i) { return data[i]; }
@ -42,8 +46,14 @@ class BaseData
explicit operator bool() const { return N() != 0; }
void SetUnit(const MString& str) { unit = str; }
void SetStandartName(const MString& str) { stname = str; }
void SetLongName(const MString& str) { lname = str; }
void SetComment(const MString& str) { comment = str; }
const MString& Unit() const { return unit; }
const MString& StandartName() const { return stname; }
const MString& LongName() const { return lname; }
const MString& Comment() const { return comment; }
};
class UngriddedData: public BaseData
@ -80,12 +90,12 @@ template<class Data> class DefaultAverager: public Data
}
public:
void Add(const Data& d)
void Add(Data&& d)
{
if(!d) return;
if(!*this)
{
*static_cast<Data*>(this) = d;
*static_cast<Data*>(this) = std::move(d);
count.resize(Data::N());
for(size_t i = 0; i < Data::N(); i++) count[i] = DataOk(this, i) ? 1 : 0;
return;

2
include/layereddata.h

@ -135,7 +135,7 @@ class LayeredData: public NCFuncs
std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const;
Data Read(const MString& vname, const BaseParameters* ip, size_t i) const;
bool Read(const MString& vname, std::map<MString, Data>& cache, const BaseParameters* ip, size_t i) const;
bool isOk() const { return nc.size() > 0; }

201
include/ncfilew.h

@ -0,0 +1,201 @@
#pragma once
#include "MString.h"
#include "traits.h"
#include <algorithm>
#include <netcdf.h>
#include <utility>
#include <vector>
using michlib::MString;
class NCFileW
{
enum Type
{
UNKNOWN,
G1V1,
G1V2,
G2V2
};
struct Var
{
MString name;
int id;
};
static constexpr auto fill = std::numeric_limits<float>::max();
int ncid;
Type type = UNKNOWN;
union
{
struct
{
int xdimid, ydimid;
};
int dimid[2];
};
int xid, yid;
int compress;
std::vector<struct Var> vars;
template<class D> static constexpr Type DetType()
{
if constexpr(ReadIs2DGeoRectArray<D>) return G1V2;
if constexpr(ReadIs2DGeoArray<D>) return G2V2;
if constexpr(ReadIs1DGeoArray<D>) return G1V1;
return UNKNOWN;
};
template<class D> static constexpr bool Is1DType() { return Is1DType(DetType<D>()); }
static constexpr bool Is1DType(Type t) { return t == G1V1; }
MString CreateFile(Type stype, const MString& name, const MString& history, int compression, size_t nx, size_t ny);
public:
static constexpr auto Fill() { return fill; }
~NCFileW() { Close(); }
template<class D> MString Create(const D& data, const MString& name, const MString& history, int compression)
{
if(type != UNKNOWN) return "File already created";
if constexpr(Is1DType<D>())
return CreateFile(DetType<D>(), name, history, compression, data.N(), 0);
else
return CreateFile(DetType<D>(), name, history, compression, data.Nx(), data.Ny());
}
size_t VarId(const MString& name) const
{
for(size_t i = 0; i < vars.size(); i++)
if(vars[i].name == name) return i + 1;
return 0;
}
bool HaveVar(const MString& name) const { return VarId(name) != 0; }
void EndDef() const
{
if(type != UNKNOWN) nc_enddef(ncid);
}
void Close()
{
if(type != UNKNOWN) nc_close(ncid);
type = UNKNOWN;
}
MString AddVariable(const MString& name, const MString& stname, const MString& lname, const MString& units, const MString& comment);
template<class D, class Op> MString WriteVariable(const D& data, size_t varid, Op op) const
{
static constexpr auto dtype = DetType<D>();
if(type == UNKNOWN) return "File not open";
if(type != dtype) return "Incompatible data type";
if(varid == 0) return "Incorrect variable";
if constexpr(dtype == UNKNOWN) return "Unknown data type";
size_t v = varid - 1;
EndDef();
int ret;
if constexpr(Is1DType(dtype))
{
const size_t i = 0, c = data.N();
float buf[c];
for(size_t ix = 0; ix < c; ix++) buf[ix] = data.IsFill(ix) ? fill : op(ix);
ret = nc_put_vara_float(ncid, vars[v].id, &i, &c, buf);
if(ret != NC_NOERR) return "Can't write " + vars[v].name + " variable in the netcdf file";
}
else
{
const size_t i[2] = {0, 0};
const size_t c[2] = {data.Nx(), data.Ny()};
float buf[c[0] * c[1]];
for(size_t ix = 0; ix < c[0]; ix++)
for(size_t iy = 0; iy < c[1]; iy++) buf[ix * c[1] + iy] = data.IsFill(ix, iy) ? fill : op(ix, iy);
ret = nc_put_vara_float(ncid, vars[v].id, i, c, buf);
if(ret != NC_NOERR) return "Can't write " + vars[v].name + " variable in the netcdf file";
}
return "";
}
template<class D, class Op> MString WriteVariable(const D& data, const MString& name, Op op) const { return WriteVariable(data, VarId(name), op); }
template<class D> MString WriteVariable(const D& data, size_t varid) const
{
if constexpr(Is1DType(DetType<D>()))
return WriteVariable(data, varid, [&data = std::as_const(data)](size_t i) { return data(i); });
else
return WriteVariable(data, varid, [&data = std::as_const(data)](size_t i, size_t j) { return data(i, j); });
}
template<class D> MString WriteVariable(const D& data, const MString& name) const { return WriteVariable(data, VarId(name)); }
template<class D> MString WriteGrid(const D& data) const
{
static constexpr auto dtype = DetType<D>();
if(type == UNKNOWN) return "File not open";
if(type != dtype) return "Incompatible data type";
EndDef();
int ret;
if constexpr(dtype == G1V1)
{
const size_t i = 0, c = data.N();
float bufx[c];
float bufy[c];
for(size_t ix = 0; ix < c; ix++)
{
bufx[ix] = data.Lon(ix);
bufy[ix] = data.Lat(ix);
}
ret = nc_put_vara_float(ncid, xid, &i, &c, bufx);
if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file";
ret = nc_put_vara_float(ncid, yid, &i, &c, bufy);
if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file";
}
else if constexpr(dtype == G1V2)
{
const size_t i = 0, cx = data.Nx(), cy = data.Ny();
float bufx[cx];
float bufy[cy];
for(size_t ix = 0; ix < cx; ix++) bufx[ix] = data.Ix2Lon(ix);
for(size_t iy = 0; iy < cy; iy++) bufy[iy] = data.Iy2Lat(iy);
ret = nc_put_vara_float(ncid, xid, &i, &cx, bufx);
if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file";
ret = nc_put_vara_float(ncid, yid, &i, &cy, bufy);
if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file";
}
else if constexpr(dtype == G2V2)
{
const size_t i[2] = {0, 0};
const size_t c[2] = {data.Nx(), data.Ny()};
float bufx[c[0] * c[1]];
float bufy[c[0] * c[1]];
for(size_t ix = 0; ix < c[0]; ix++)
for(size_t iy = 0; iy < c[1]; iy++)
{
bufx[ix * c[1] + iy] = data.Lon(ix, iy);
bufy[ix * c[1] + iy] = data.Lat(ix, iy);
}
ret = nc_put_vara_float(ncid, xid, i, c, bufx);
if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file";
ret = nc_put_vara_float(ncid, yid, i, c, bufy);
if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file";
}
else
return "Unknown data type";
return "";
}
};

35
include/simple2ddata.h

@ -18,6 +18,8 @@ class Simple2DData: public BaseData
{
}
Simple2DData CopyGrid() const { return {nx, ny, x0, y0, xstep, ystep}; }
const real& V(size_t ix, size_t iy) const { return V(iy * nx + ix); }
real& V(size_t ix, size_t iy) { return V(iy * nx + ix); }
@ -41,3 +43,36 @@ class Simple2DData: public BaseData
real XStep() const { return xstep; }
real YStep() const { return ystep; }
};
class Rect2DData: public BaseData
{
std::vector<real> lon, lat;
public:
using BaseData::IsFill;
using BaseData::V;
using BaseData::operator();
Rect2DData() = default;
Rect2DData(const std::vector<real>& lons, const std::vector<real>& lats, MString&& unit = ""): BaseData(lons.size() * lats.size(), std::move(unit)), lon(lons), lat(lats) {}
const real& V(size_t ix, size_t iy) const { return V(iy * Nx() + ix); }
real& V(size_t ix, size_t iy) { return V(iy * Nx() + ix); }
const real& operator()(size_t ix, size_t iy) const { return V(iy * Nx() + ix); }
real& operator()(size_t ix, size_t iy) { return V(iy * Nx() + ix); }
size_t Nx() const { return lon.size(); }
size_t Ny() const { return lat.size(); }
bool IsFill(size_t ix, size_t iy) const { return IsFill(iy * Nx() + ix); }
real Lon(size_t ix, [[maybe_unused]] size_t iy) const { return lon[ix]; }
real Lat([[maybe_unused]] size_t ix, size_t iy) const { return lat[iy]; }
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 Ix2Lon(size_t ix) const { return lon[ix]; }
real Iy2Lat(size_t iy) const { return lat[iy]; }
};

61
include/traits.h

@ -3,6 +3,7 @@
#include <concepts>
class BaseParameters;
using michlib::real;
template<class T>
@ -12,6 +13,13 @@ concept InfoSupported = requires(T t) {
} -> std::convertible_to<MString>;
};
template<class T>
concept HasDefVars = requires(T t) {
{
t.DefaultVars()
} -> std::convertible_to<MString>;
};
template<class T>
concept ParametersRequiredRegion = requires(T t, michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& r) {
{
@ -30,16 +38,16 @@ template<class T>
concept ParametersSupported = ParametersRequiredRegion<T> || ParametersNotRequiredRegion<T>;
template<class T>
concept ReadPSupported = requires(T t, const MString& vname, const BaseParameters* ip, size_t i) {
concept ReadPSupported = requires(T t, const MString& vname, std::map<MString, typename T::Data>& cache, const BaseParameters* ip, size_t i) {
{
t.Read(vname, ip, i)(0)
} -> std::convertible_to<real>;
t.Read(vname, cache, ip, i)
} -> std::same_as<bool>;
};
template<class T>
concept ReadSupported = requires(T t, const MString& vname, size_t i) {
concept ReadSupported = requires(T t, const MString& vname, std::map<MString, typename T::Data>& cache, size_t i) {
{
t.Read(vname, i)(0)
} -> std::convertible_to<real>;
t.Read(vname, cache, i)
} -> std::same_as<bool>;
};
template<class T>
@ -69,17 +77,16 @@ concept HaveXYStep = requires(T t) {
template<typename...> using void_ = void;
template<class D, bool RP, bool R> struct GetReadType_;
template<class D, bool R> struct GetReadType_;
template<class D, bool R> struct GetReadType_<D, true, R>
template<class D> struct GetReadType_<D, true>
{
using p = BaseParameters*;
using type = decltype(D().Read(MString(), p(), 0));
using type = typename D::Data;
};
template<class D> struct GetReadType_<D, false, true>
template<class D> struct GetReadType_<D, false>
{
using type = decltype(D().Read(MString(), 0));
using type = std::decay_t<D>;
};
template<class T>
@ -119,24 +126,44 @@ template<class Act, class S> consteval bool IsDisabled()
}
}
template<class D> using ReadType = internal::GetReadType_<D, ReadPSupported<D>, ReadSupported<D>>::type;
template<class D> using ReadType = internal::GetReadType_<D, ReadPSupported<D> || ReadSupported<D>>::type;
template<class T>
concept ReadIsGrid = requires {
{
ReadType<T>().XStep()
std::declval<ReadType<T>>().XStep()
} -> std::convertible_to<real>;
{
ReadType<T>().YStep()
std::declval<ReadType<T>>().YStep()
} -> std::convertible_to<real>;
};
template<class T>
concept ReadIs2DGeoRectArray = requires {
{
ReadType<T>().Ix2Lon(0)
std::declval<ReadType<T>>().Ix2Lon(0)
} -> std::convertible_to<real>;
{
std::declval<ReadType<T>>().Iy2Lat(0)
} -> std::convertible_to<real>;
};
template<class T>
concept ReadIs2DGeoArray = requires {
{
std::declval<ReadType<T>>().Lon(0, 0)
} -> std::convertible_to<real>;
{
std::declval<ReadType<T>>().Lat(0, 0)
} -> std::convertible_to<real>;
};
template<class T>
concept ReadIs1DGeoArray = requires {
{
std::declval<ReadType<T>>().Lon(0)
} -> std::convertible_to<real>;
{
ReadType<T>().Iy2Lat(0)
std::declval<ReadType<T>>().Lat(0)
} -> std::convertible_to<real>;
};

22
include/uvdata.h

@ -21,6 +21,7 @@ class UVDataStorage
std::vector<real> u, v, u2;
size_t nx = 0, ny = 0;
Metric metric;
MString unit;
UVDataStorage() = default;
real D(real lon1, real lat1, real lon2, real lat2) const
@ -38,6 +39,8 @@ class UVDataStorage
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); }
void SetUnit(const MString& str) { unit = str; }
public:
enum STPOINT
{
@ -82,7 +85,22 @@ class UVDataStorage
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; }
bool IsFill(size_t i) const { return IsCoast(i); }
bool IsFill(size_t ix, size_t iy) const { return IsCoast(ix, iy); }
static real Fillval() { return fillval; }
const MString& Unit() const { return unit; }
const MString DUnit() const
{
switch(metric)
{
case(Metric::EUCLID): return "";
case(Metric::SPHERIC): return "km";
}
return "";
}
};
template<class OneVarData, bool isgrid> class UVDataDims;
@ -128,6 +146,7 @@ template<class OneVarData> class UVDataDims<OneVarData, true>: public UVDataStor
U2(i) = u(i) * u(i) + v(i) * v(i);
}
}
SetUnit(u.Unit());
}
public:
@ -138,6 +157,9 @@ template<class OneVarData> class UVDataDims<OneVarData, true>: public UVDataStor
real XStep() const { return xstep; }
real YStep() const { return ystep; }
real Ix2Lon(size_t ix) const { return x0 + ix * xstep; }
real Iy2Lat(size_t iy) const { return y0 + iy * ystep; }
};
template<class OneVarData> class UVDataDims<OneVarData, false>: public UVDataStorage

43
sources/AVISOLOCAL.cpp

@ -86,9 +86,10 @@ std::pair<const BaseParameters*, MString> AVISOLOCALData::Parameters(michlib_int
return {ppar.release(), ""};
}
AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParameters* ip, size_t i) const
bool AVISOLOCALData::Read(const MString& vname, std::map<MString, AVISOLOCALData::Data>& cache, const BaseParameters* ip, size_t i) const
{
if(!isOk()) return Data();
if(cache.contains(vname)) return true;
if(!isOk()) return false;
auto p = dynamic_cast<const struct Parameters*>(ip);
NCFileA nc;
@ -96,7 +97,7 @@ AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParame
bool isfloat = false, isint2 = false, isint4 = false;
nc.Reset(datapath + "/uv-" + times[i].ToString() + ".nc");
if(!nc) return Data();
if(!nc) return false;
{
auto head = nc.Header();
for(const auto& v: head.Variables())
@ -119,28 +120,36 @@ AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParame
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;
out.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
for(size_t ind = 0; ind < out.N(); ind++)
if(!(Read("u", cache, ip, i) && Read("v", cache, ip, i))) return false;
cache[vname] = cache.at("u").CopyGrid();
auto& U = cache.at(vname);
const auto& u = cache.at("u");
const auto& v = cache.at("v");
if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
for(size_t ind = 0; ind < U.N(); ind++)
{
if(u.IsFill(ind) || v.IsFill(ind))
out.V(ind) = out.Fillval();
U.V(ind) = U.Fillval();
else
out.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
}
return out;
return true;
}
return Data();
return false;
}
// Direct read
if(isint2) return ReadVarRaw<int2>(nc, name, p);
if(isint4) return ReadVarRaw<int4>(nc, name, p);
if(isfloat) return ReadVarRaw<float>(nc, name, p);
return Data();
Data data;
if(isint2) data = ReadVarRaw<int2>(nc, name, p);
if(isint4) data = ReadVarRaw<int4>(nc, name, p);
if(isfloat) data = ReadVarRaw<float>(nc, name, p);
if(data)
{
cache[vname] = std::move(data);
return true;
}
return false;
}
template<class DataType> AVISOLOCALData::Data AVISOLOCALData::ReadVarRaw(const NCFileA& nc, const MString& name, const struct AVISOLOCALData::Parameters* p) const

2
sources/AVISOLOCAL.h

@ -60,7 +60,7 @@ class AVISOLOCALData: public NCFuncs
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;
bool Read(const MString& vname, std::map<MString, Data>& cache, const BaseParameters* ip, size_t i) const;
template<class DataType> Data ReadVarRaw(const NCFileA& nc, const MString& name, const struct Parameters* p) const;
};

32
sources/BINFILE.cpp

@ -53,9 +53,10 @@ MString BINFILEData::Open(const CLArgs& args)
return "";
}
BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const
bool BINFILEData::Read(const MString& vname, std::map<MString, BINFILEData::Data>& cache, size_t i) const
{
if(!isOk()) return Data();
if(cache.contains(vname)) return true;
if(!isOk()) return false;
// Only rectangular grids are supported
real xs = (data->Lon(data->Nx() - 1, data->Ny() - 1) - data->Lon(0, 0)) / (data->Nx() - 1);
@ -66,19 +67,20 @@ BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const
if(vname == "U" || vname == "U2")
{
bool square = vname == "U2";
auto u = Read("u", i);
auto v = Read("v", i);
if(!(u && v)) return Data();
auto out = u;
out.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
for(size_t ind = 0; ind < out.N(); ind++)
if(!(Read("u", cache, i) && Read("v", cache, i))) return false;
cache[vname] = cache.at("u").CopyGrid();
auto& U = cache.at(vname);
const auto& u = cache.at("u");
const auto& v = cache.at("v");
if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
for(size_t ind = 0; ind < U.N(); ind++)
{
if(u.IsFill(ind) || v.IsFill(ind))
out.V(ind) = out.Fillval();
U.V(ind) = U.Fillval();
else
out.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
}
return out;
return true;
}
if(vname == "u" || vname == "v")
@ -92,7 +94,13 @@ BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const
else
out.V(ix, iy) = out.Fillval();
}
if(out)
{
cache[vname] = std::move(out);
return true;
}
}
return out;
return false;
}

2
sources/BINFILE.h

@ -44,5 +44,5 @@ class BINFILEData
bool CheckVar(const MString& vname) const { return vname == "u" || vname == "v" || vname == "U" || vname == "U2"; }
Data Read(const MString& vname, size_t i) const;
bool Read(const MString& vname, std::map<MString, Data>& cache, size_t i) const;
};

12
sources/MODISBINLOCAL.cpp

@ -83,8 +83,9 @@ std::pair<const BaseParameters*, MString> MODISBINLOCALData::Parameters(michlib_
return {ppar.release(), ""};
}
MODISBINLOCALData::Data MODISBINLOCALData::Read(const MString& var, const BaseParameters* ip, size_t tind) const
bool MODISBINLOCALData::Read(const MString& vname, std::map<MString, MODISBINLOCALData::Data>& cache, const BaseParameters* ip, size_t tind) const
{
if(cache.contains(vname)) return true;
auto p = pointer_cast<const struct Parameters*>(ip);
Averager<Data> out;
@ -92,9 +93,14 @@ MODISBINLOCALData::Data MODISBINLOCALData::Read(const MString& var, const BasePa
for(const auto& suf: dataset)
{
Data dat = ReadFile(suf, p, tind);
if(dat) out.Add(dat);
if(dat) out.Add(std::move(dat));
}
return out.Div();
if(out)
{
cache[vname] = std::move(out.Div());
return true;
}
return false;
}
MODISBINLOCALData::Data MODISBINLOCALData::ReadFile(const MString& suf, const struct MODISBINLOCALData::Parameters* p, size_t tind) const

11
sources/MODISBINLOCAL.h

@ -40,9 +40,18 @@ class MODISBINLOCALData
MString Info() const;
MString Open(const CLArgs& args);
Data Read(const MString& var, const BaseParameters* ip, size_t tind) const;
bool Read(const MString& vname, std::map<MString, Data>& cache, const BaseParameters* ip, size_t tind) const;
Data ReadFile(const MString& suf, const struct Parameters* p, size_t tind) const;
MString DefaultVars() const
{
if(dataset.size() == 0) return "";
if(dataset[0] == "CHL" || dataset[0] == "T_CHL" || dataset[0] == "A_CHL")
return "chl";
else
return "temp";
}
std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const;
bool CheckVar(const MString& vname) const

526
sources/VYLET.cpp

@ -0,0 +1,526 @@
#define MICHLIB_NOSOURCE
#include "VYLET.h"
MString VYLETData::Info() const
{
if(!vylet) return "";
MString out;
michlib::CompiledParser xy2lon, xy2lat;
real x, y;
michlib::ParserVars pv;
pv["x"] = &x;
pv["y"] = &y;
vylet->UsePrefix("Datafile_Info");
MString xy2lonstr = vylet->ParameterSValue("xy2lon", "%y");
MString xy2latstr = vylet->ParameterSValue("xy2lat", "%x");
ArifmeticCompiler(xy2lonstr, xy2lon, &pv);
ArifmeticCompiler(xy2latstr, xy2lat, &pv);
real lonb, latb, lone, late;
vylet->UsePrefix("");
x = vylet->ParameterRValue("x0", 0.0);
y = vylet->ParameterRValue("y0", 0.0);
xy2lon.Run(lonb);
xy2lat.Run(latb);
x = vylet->ParameterRValue("x1", 0.0);
y = vylet->ParameterRValue("y1", 0.0);
xy2lon.Run(lone);
xy2lat.Run(late);
out += "Start time: " + start.ToString() + "\n";
out += "End time: " + end.ToString() + " " + (invtime ? "(backward integration)" : "(forward integration)") + "\n";
out += "Region: (" + MString(lonb) + " : " + lone + ") x (" + latb + " : " + late + ")\n";
out += "Grid:" + MString(" ") + lons.size() + "x" + lats.size() + "\n";
out += HasCoast() ? "Coast: checked\n" : "Coast: not checked\n";
out += LengthFast() ? "Trajectory length: fast calculation\n" : "Trajectory length: exact calculation\n";
out += "Borders: ";
bool needcomma = false;
if(Left() >= 0.0)
{
real lon;
if(needcomma) out += ", ";
out += "left=";
x = Left();
xy2lon.Run(lon);
out += lon;
needcomma = true;
}
if(Right() <= maxx)
{
real lon;
if(needcomma) out += ", ";
out += "right=";
x = Right();
xy2lon.Run(lon);
out += lon;
needcomma = true;
}
if(Down() >= 0.0)
{
real lat;
if(needcomma) out += ", ";
out += "bottom=";
y = Down();
xy2lat.Run(lat);
out += lat;
needcomma = true;
}
if(Up() <= maxy)
{
real lat;
if(needcomma) out += ", ";
out += "bottom=";
y = Up();
xy2lat.Run(lat);
out += lat;
needcomma = true;
}
if(!needcomma) out += "no";
out += "\n";
MString vars = "vylD";
if(HasLyap()) vars += ", vylL";
if(HasTime()) vars += ", vylT";
vars += ", vylNx, vylNy, vylRx, vylRy, vylEx, vylEy, vylPhip, vylPhim, vylPhit";
if(HasLyap()) vars += ", vyldS";
vars += ", vylAngle, vylLen";
if(HasTime()) vars += ", vylTmask";
out += "Supported variables: " + vars + "\n";
return out;
}
MString VYLETData::Open(const CLArgs& args)
{
if(!args.contains("dataset")) return "path to data not specified";
MString dataset = args.at("dataset");
decltype(vylet) nvylet;
michlib::RegExpSimple havex("%x"), havey("%y");
nvylet.reset(new michlib::BFileR);
if(nvylet->Open(dataset) != ERR_NOERR) return "Can't open file " + dataset;
nvylet->UsePrefix("ProgramInfo");
if(nvylet->ParameterSValue("Task", "") != "AdvInt:Vylet") return "File " + dataset + " is not vylet file";
nvylet->UsePrefix("");
if(nvylet->ParameterSValue("nettype", "") != "SQUARE") return "File " + dataset + " have unsupported net type";
MString method = nvylet->ParameterSValue("Method", "");
if(!(method == "Bicubic" || method == "BicubicL" || method == "BicubicI" || method == "BicubicIL")) return "File " + dataset + " have unsupported integration method";
invtime = (method == "BicubicI" || method == "BicubicIL");
nvylet->UsePrefix("Datafile_Info");
MString xy2lonstr = nvylet->ParameterSValue("xy2lon", "%y");
MString xy2latstr = nvylet->ParameterSValue("xy2lat", "%x");
if(havey.Match(xy2lonstr.Buf()) || havex.Match(xy2latstr.Buf())) return "File " + dataset + " have unsupported grid";
auto res = ref.FromString(nvylet->ParameterSValue(invtime ? "EndDate" : "BeginDate", ""));
if(!res) return "Can't read reference time";
nvylet->UsePrefix("");
start = R2Time(nvylet->ParameterRValue("tbeg", 0.0));
end = R2Time(nvylet->ParameterRValue("tmax", 0.0));
nvylet->UsePrefix("Datafile");
auto dx = nvylet->ParameterRValue("dx", 0.0);
auto dy = nvylet->ParameterRValue("dy", 0.0);
auto nx = nvylet->ParameterUValue("nx", 0);
auto ny = nvylet->ParameterUValue("ny", 0);
maxx = dx * (nx - 1);
maxy = dy * (ny - 1);
michlib::CompiledParser xy2lon, xy2lat;
real x, y;
michlib::ParserVars pv;
pv["x"] = &x;
pv["y"] = &y;
ArifmeticCompiler(xy2lonstr, xy2lon, &pv);
ArifmeticCompiler(xy2latstr, xy2lat, &pv);
nvylet->UsePrefix("");
auto nlon = nvylet->ParameterUValue("Nx", 0) + 1;
auto nlat = nvylet->ParameterUValue("Ny", 0) + 1;
lons.resize(nlon);
lats.resize(nlat);
elon.resize(nlon * nlat);
elat.resize(nlon * nlat);
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++)
{
x = (*nvylet)[2][iy * lons.size() + ix];
y = (*nvylet)[3][iy * lons.size() + ix];
xy2lon.Run(elon[iy * lons.size() + ix]);
xy2lat.Run(elat[iy * lons.size() + ix]);
}
for(size_t ix = 0; ix < nlon; ix++)
{
x = (*nvylet)[0][ix];
y = (*nvylet)[1][ix];
xy2lon.Run(lons[ix]);
}
for(size_t iy = 0; iy < nlat; iy++)
{
x = (*nvylet)[0][iy * nlon];
y = (*nvylet)[1][iy * nlon];
xy2lat.Run(lats[iy]);
}
vylet = std::move(nvylet);
return "";
}
bool VYLETData::Read(const MString& vname, std::map<MString, VYLETData::Data>& cache, size_t tind) const
{
if(tind != 0) return false;
if(cache.contains(vname)) return true;
Data out;
if(vname == "vylD") out = ReadD();
if(vname == "vylL") out = ReadL();
if(vname == "vylT") out = ReadT();
if(vname == "vylNx") out = ReadNx();
if(vname == "vylNy") out = ReadNy();
if(vname == "vylRx") out = ReadRx();
if(vname == "vylRy") out = ReadRy();
if(vname == "vylEx") out = ReadEx();
if(vname == "vylEy") out = ReadEy();
if(vname == "vylPhip") out = ReadPhip();
if(vname == "vylPhim") out = ReadPhim();
if(vname == "vylPhit") out = ReadPhit();
if(vname == "vyldS") out = ReaddS();
if(vname == "vylAngle") out = ReadAngle();
if(vname == "vylLen") out = ReadLen();
if(vname == "vylTmask") out = ReadTmask();
if(!out) return false;
cache[vname] = std::move(out);
return true;
}
VYLETData::Data VYLETData::ReadD() const
{
Data out(lons, lats);
if(!vylet) return Data();
const real xl = Left();
const real xr = Right();
const real yd = Down();
const real yu = Up();
real x, y;
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++)
{
out(ix, iy) = 0.0;
x = (*vylet)[2][iy * lons.size() + ix];
y = (*vylet)[3][iy * lons.size() + ix];
if(yd > 0.0 && y < yd) out(ix, iy) = -1.0;
if(xl > 0.0 && x < xl) out(ix, iy) = -2.0;
if(yu < maxy && y > yu) out(ix, iy) = -3.0;
if(xr < maxx && x > xr) out(ix, iy) = -4.0;
if(out(ix, iy) >= 0.0)
out(ix, iy) = 6371.0 * michlib::GCD(M_PI * lons[ix] / 180.0, M_PI * lats[iy] / 180.0, M_PI * elon[iy * lons.size() + ix] / 180.0, M_PI * elat[iy * lons.size() + ix] / 180.0);
}
out.SetUnit("km");
out.SetLongName("Distance between start and end points");
out.SetComment("Special values: -1.0 - trajectory leave region via south border, -2.0 - trajectory leave region via west border, -3.0 - trajectory leave region via north border, "
"-4.0 - trajectory leave region via east border");
return out;
}
VYLETData::Data VYLETData::ReadL() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto lcol = Name2ColNum("Log(G.sigma1)");
auto tcol = Name2ColNum("time");
if(lcol == 0 || tcol == 0) return Data();
MDateTime time;
real lambda, days;
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++)
{
time = R2Time((*vylet)[tcol - 1][iy * lons.size() + ix]);
lambda = (*vylet)[lcol - 1][iy * lons.size() + ix];
days = static_cast<real>(time - start) / MDateTime::secondsperday;
out(ix, iy) = lambda / days;
}
out.SetUnit("days-1");
out.SetLongName("Lyapunov exponent");
return out;
}
VYLETData::Data VYLETData::ReadT() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto tcol = Name2ColNum("time");
if(tcol == 0) return Data();
const real xl = Left();
const real xr = Right();
const real yd = Down();
const real yu = Up();
vylet->UsePrefix("");
const real acc = vylet->ParameterRValue("accuracy", 1.0);
const real tstep = 2.0 * M_PI * 1000.0 / acc;
const real maxdays = static_cast<real>(end - start) / MDateTime::secondsperday;
MDateTime time;
real days;
real x, y;
bool inside, maxtime;
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++)
{
x = (*vylet)[2][iy * lons.size() + ix];
y = (*vylet)[3][iy * lons.size() + ix];
time = R2Time((*vylet)[tcol - 1][iy * lons.size() + ix]);
days = static_cast<real>(time - start) / MDateTime::secondsperday;
if(days <= tstep * 1.5) days = 0.0;
inside = x > xl && x < xr && y > yd && y < yu;
maxtime = days >= maxdays - tstep * 0.5;
if(maxtime) days = maxdays;
if(inside && !maxtime) days = -days;
out(ix, iy) = days;
}
out.SetUnit("days");
out.SetLongName("Time to reach border or coast");
out.SetComment("Special values: 0.0 - initial point on the land, " + MString(maxdays) + " - trajectory don't reach border, negative values - trajectory beached on the coast");
return out;
}
VYLETData::Data VYLETData::ReadNx() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto ncol = Name2ColNum("nx=0");
if(ncol == 0) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = (*vylet)[ncol - 1][iy * lons.size() + ix];
out.SetLongName("Number of moments u=0");
return out;
}
VYLETData::Data VYLETData::ReadNy() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto ncol = Name2ColNum("ny=0");
if(ncol == 0) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = (*vylet)[ncol - 1][iy * lons.size() + ix];
out.SetLongName("Number of moments v=0");
return out;
}
VYLETData::Data VYLETData::ReadRx() const
{
Data out(lons, lats);
if(!vylet) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = elon[iy * lons.size() + ix] - lons[ix];
out.SetLongName("Displacement in zonal direction");
return out;
}
VYLETData::Data VYLETData::ReadRy() const
{
Data out(lons, lats);
if(!vylet) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = elat[iy * lons.size() + ix] - lats[iy];
out.SetLongName("Displacement in meridional direction");
return out;
}
VYLETData::Data VYLETData::ReadEx() const
{
Data out(lons, lats);
if(!vylet) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = elon[iy * lons.size() + ix];
out.SetLongName("Final longitude");
return out;
}
VYLETData::Data VYLETData::ReadEy() const
{
Data out(lons, lats);
if(!vylet) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = elat[iy * lons.size() + ix];
out.SetLongName("Final latitude");
return out;
}
VYLETData::Data VYLETData::ReadPhip() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto ncol = Name2ColNum("nphip");
if(ncol == 0) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = (*vylet)[ncol - 1][iy * lons.size() + ix];
out.SetLongName("Number of counterclockwise rotations");
return out;
}
VYLETData::Data VYLETData::ReadPhim() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto ncol = Name2ColNum("nphim");
if(ncol == 0) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = (*vylet)[ncol - 1][iy * lons.size() + ix];
out.SetLongName("Number of clockwise rotations");
return out;
}
VYLETData::Data VYLETData::ReadPhit() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto pcol = Name2ColNum("nphip");
auto mcol = Name2ColNum("nphim");
if(pcol == 0 || mcol == 0) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = (*vylet)[pcol - 1][iy * lons.size() + ix] - (*vylet)[mcol - 1][iy * lons.size() + ix];
out.SetLongName("Difference between the number of rotations counterclockwise and clockwise");
return out;
}
VYLETData::Data VYLETData::ReaddS() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto l1col = Name2ColNum("Log(G.sigma1)");
auto l2col = Name2ColNum("Log(G.sigma2)");
if(l1col == 0 || l2col == 0) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = (*vylet)[l1col - 1][iy * lons.size() + ix] + (*vylet)[l2col - 1][iy * lons.size() + ix];
out.SetLongName("Logarithm of area change multiplier");
return out;
}
VYLETData::Data VYLETData::ReadAngle() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto ncol = Name2ColNum("angle");
if(ncol == 0) return Data();
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = (*vylet)[ncol - 1][iy * lons.size() + ix] / (2.0 * M_PI);
out.SetLongName("Total rotation angle normalized to 2π");
return out;
}
VYLETData::Data VYLETData::ReadLen() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto ncol = Name2ColNum("length");
if(ncol == 0) return Data();
bool sisangle = !LengthFast();
real mul = sisangle ? (6371.0 * M_PI / (60.0 * 180.0)) : 1.0;
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++) out(ix, iy) = mul * (*vylet)[ncol - 1][iy * lons.size() + ix];
if(sisangle)
{
out.SetUnit("km");
out.SetLongName("Trajectory length");
}
out.SetLongName("Trajectory length (in abstract units)");
return out;
}
VYLETData::Data VYLETData::ReadTmask() const
{
Data out(lons, lats);
if(!vylet) return Data();
auto tcol = Name2ColNum("time");
if(tcol == 0) return Data();
vylet->UsePrefix("");
const real acc = vylet->ParameterRValue("accuracy", 1.0);
const real tstep = 2.0 * M_PI * 1000.0 / acc;
const real maxdays = static_cast<real>(end - start) / MDateTime::secondsperday;
MDateTime time;
real days;
bool maxtime;
for(size_t iy = 0; iy < lats.size(); iy++)
for(size_t ix = 0; ix < lons.size(); ix++)
{
time = R2Time((*vylet)[tcol - 1][iy * lons.size() + ix]);
days = static_cast<real>(time - start) / MDateTime::secondsperday;
maxtime = days >= maxdays - tstep * 0.5;
out(ix, iy) = maxtime ? 1.0 : NAN;
}
out.SetLongName("Flag");
out.SetComment("Values: 1.0 - the trajectory did not reach either the coast or the borders, NaN - overwise");
return out;
}

145
sources/VYLET.h

@ -0,0 +1,145 @@
#pragma once
#include "BFileR.h"
#include "ParseArgs.h"
#include "mdatetime.h"
#include "mregex.h"
#include "simple2ddata.h"
#include <memory>
#include <set>
using michlib::M_PI;
using michlib::MDateTime;
using michlib::real;
using michlib::Round;
class VYLETData
{
std::unique_ptr<michlib::BFileR> vylet;
MDateTime start, end, ref;
bool invtime;
real maxx, maxy;
std::vector<real> lons, lats;
std::vector<real> elon, elat;
MDateTime R2Time(real r) const
{
auto sec = static_cast<time_t>(Round(r * MDateTime::secondsperday));
return ref + (invtime ? -sec : sec);
};
public:
using Data = Rect2DData;
private:
auto Name2ColNum(const MString& name) const
{
auto nc = vylet->Columns();
decltype(nc) i;
for(i = 1; i <= nc; i++)
if(vylet->ColumnName(i) == name) return i;
i = 0;
return i;
}
Data ReadD() const;
Data ReadL() const;
Data ReadT() const;
Data ReadNx() const;
Data ReadNy() const;
Data ReadRx() const;
Data ReadRy() const;
Data ReadEx() const;
Data ReadEy() const;
Data ReadPhip() const;
Data ReadPhim() const;
Data ReadPhit() const;
Data ReaddS() const;
Data ReadAngle() const;
Data ReadLen() const;
Data ReadTmask() const;
public:
static constexpr const char* name = "VYLET";
static constexpr const char* disabledactions = "genintfile uv";
MString DefaultVars() const
{
MString vars = "vylD,vylRx,vylRy,vylEx,vylEy,vylAngle,vylLen";
if(HasLyap()) vars += ",vylL";
if(HasTime()) vars += ",vylT";
return vars;
}
bool CheckVar(const MString& vname) const
{
std::set<MString> vars{"vylD", "vylNx", "vylNy", "vylRx", "vylRy", "vylEx", "vylEy", "vylPhip", "vylPhim", "vylPhit", "vylAngle", "vylLen"};
if(HasLyap()) vars.insert("vylL");
if(HasLyap()) vars.insert("vyldS");
if(HasTime()) vars.insert("vylT");
if(HasTime()) vars.insert("vylTmask");
return vars.contains(vname);
}
bool Read(const MString& vname, std::map<MString, Data>& cache, size_t tind) const;
size_t NTimes() const { return vylet ? 1 : 0; }
MDateTime Time(size_t i) const
{
if(i == 0 && vylet) return start;
return MDateTime();
}
bool HasLyap() const
{
if(!vylet) return false;
vylet->UsePrefix("");
MString method = vylet->ParameterSValue("Method", "");
return method == "BicubicL" || method == "BicubicIL";
}
bool HasCoast() const
{
if(!vylet) return false;
vylet->UsePrefix("");
return vylet->ParameterBValue("checkcoast", true);
}
real Left() const
{
vylet->UsePrefix("");
return vylet->ParameterRValue("xl", -1.0);
}
real Right() const
{
vylet->UsePrefix("");
return vylet->ParameterRValue("xr", maxx + 1.0);
}
real Down() const
{
vylet->UsePrefix("");
return vylet->ParameterRValue("yd", -1.0);
}
real Up() const
{
vylet->UsePrefix("");
return vylet->ParameterRValue("yu", maxy + 1.0);
}
bool HasBorders() const { return Left() >= 0.0 || Right() <= maxx || Down() >= 0.0 || Up() <= maxy; }
bool HasTime() const { return HasCoast() || HasBorders(); }
bool LengthFast() const
{
vylet->UsePrefix("");
return !vylet->ParameterBValue("sisangle", false);
}
MString Info() const;
MString Open(const CLArgs& args);
};

7
src/ParseArgs.cpp

@ -33,3 +33,10 @@ CLArgs ParseArgs(int argc, char** argv)
return out;
}
MString GetExt(const MString& fname)
{
for(size_t i = fname.Len(); i != 0; i--)
if(fname[i - 1] == '.') return fname.SubStr(i + 1, fname.Len() - i);
return "";
}

105
src/layereddata.cpp

@ -178,10 +178,10 @@ std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_intern
return {ppar.release(), ""};
}
LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters* ip, size_t i) const
bool LayeredData::Read(const MString& vname, std::map<MString, LayeredData::Data>& cache, const BaseParameters* ip, size_t i) const
{
if(!isOk()) return Data();
bool nodepth = false;
if(cache.contains(vname)) return true;
if(!isOk()) return false;
auto p = dynamic_cast<const struct Parameters*>(ip);
auto [name, id, tid] = VarNameLoc(vname, times[i]);
@ -190,90 +190,101 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters*
// ptemp from temp and sal
if(vname == "ptemp")
{
auto temp = Read("temp", ip, i);
auto sal = Read("sal", ip, i);
if(!(temp && sal)) return Data();
auto out = temp;
out.SetUnit("degrees_C");
for(size_t ind = 0; ind < out.N(); ind++)
if(!(Read("temp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
cache[vname] = cache.at("temp").CopyGrid();
auto& ptemp = cache.at(vname);
const auto& temp = cache.at("temp");
const auto& sal = cache.at("sal");
ptemp.SetUnit("degrees_C");
for(size_t ind = 0; ind < ptemp.N(); ind++)
{
if(temp.IsFill(ind) || sal.IsFill(ind))
out.V(ind) = out.Fillval();
ptemp.V(ind) = ptemp.Fillval();
else
out.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
ptemp.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), Depth(p->layer), ptemp.Lon(ind), ptemp.Lat(ind));
}
return out;
return true;
}
// temp from ptemp and sal
if(vname == "temp")
{
auto temp = Read("ptemp", ip, i);
auto sal = Read("sal", ip, i);
if(!(temp && sal)) return Data();
auto out = temp;
out.SetUnit("degrees_C");
for(size_t ind = 0; ind < out.N(); ind++)
if(!(Read("ptemp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
cache[vname] = cache.at("ptemp").CopyGrid();
auto& temp = cache.at(vname);
const auto& ptemp = cache.at("ptemp");
const auto& sal = cache.at("sal");
temp.SetUnit("degrees_C");
for(size_t ind = 0; ind < temp.N(); ind++)
{
if(temp.IsFill(ind) || sal.IsFill(ind))
out.V(ind) = out.Fillval();
if(ptemp.IsFill(ind) || sal.IsFill(ind))
temp.V(ind) = temp.Fillval();
else
out.V(ind) = PTemp2Temp(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
temp.V(ind) = PTemp2Temp(ptemp.V(ind), sal.V(ind), Depth(p->layer), temp.Lon(ind), temp.Lat(ind));
}
return out;
return true;
}
// pdens from temp and sal
if(vname == "pdens")
{
bool tempispot = HaveVar("ptemp");
auto temp = Read(tempispot ? "ptemp" : "temp", ip, i);
auto sal = Read("sal", ip, i);
if(!(temp && sal)) return Data();
auto out = temp;
out.SetUnit("kg m-3");
for(size_t ind = 0; ind < out.N(); ind++)
if(!(Read(tempispot ? "ptemp" : "temp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
cache[vname] = cache.at("sal").CopyGrid();
auto& pdens = cache.at(vname);
const auto& temp = cache.at(tempispot ? "ptemp" : "temp");
const auto& sal = cache.at("sal");
pdens.SetUnit("kg m-3");
for(size_t ind = 0; ind < pdens.N(); ind++)
{
if(temp.IsFill(ind) || sal.IsFill(ind))
out.V(ind) = out.Fillval();
pdens.V(ind) = pdens.Fillval();
else
out.V(ind) = tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind))
: Temp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
pdens.V(ind) = tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), pdens.Lon(ind), pdens.Lat(ind))
: Temp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), pdens.Lon(ind), pdens.Lat(ind));
}
return out;
return true;
}
// 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;
out.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
for(size_t ind = 0; ind < out.N(); ind++)
if(!(Read("u", cache, ip, i) && Read("v", cache, ip, i))) return false;
cache[vname] = cache.at("u").CopyGrid();
auto& U = cache.at(vname);
const auto& u = cache.at("u");
const auto& v = cache.at("v");
if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
for(size_t ind = 0; ind < U.N(); ind++)
{
if(u.IsFill(ind) || v.IsFill(ind))
out.V(ind) = out.Fillval();
U.V(ind) = U.Fillval();
else
out.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
}
return out;
return true;
}
return Data();
return false;
}
// Direct read
bool nodepth = false;
Data data;
auto head = nc[id]->Header();
for(const auto& v: head.Variables())
if(v.Name() == name)
{
if(v.Dimensions().size() == 3) nodepth = true;
if(v.Type().Id() == NC_SHORT) return ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_INT) return ReadVarRaw<int>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_DOUBLE) return ReadVarRaw<double>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_SHORT) data = ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_INT) data = ReadVarRaw<int>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_FLOAT) data = ReadVarRaw<float>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_DOUBLE) data = ReadVarRaw<double>(nc[id], name, tid, nodepth, p);
if(data)
{
cache[vname] = std::move(data);
return true;
}
}
return Data();
return false;
}
template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct LayeredData::Parameters* p) const

105
src/ncfilew.cpp

@ -0,0 +1,105 @@
#define MICHLIB_NOSOURCE
#include "ncfilew.h"
MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, const MString& history, int compression, size_t nx, size_t ny)
{
if(stype == UNKNOWN) return "Can't determine file type";
compress = compression;
int ret;
MString text;
ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &ncid);
if(ret != NC_NOERR) return "Can't create netcdf file: " + name;
ret = nc_put_att_text(ncid, NC_GLOBAL, "history", history.Len() + 1, history.Buf());
if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
switch(stype)
{
case(G1V1):
{
ret = nc_def_dim(ncid, "i", nx, &xdimid);
if(ret != NC_NOERR) return "Can't create dimension in the netcdf file";
ret = nc_def_var(ncid, "longitude", NC_FLOAT, 1, &xdimid, &xid);
if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file";
ret = nc_def_var(ncid, "latitude", NC_FLOAT, 1, &xdimid, &yid);
if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file";
break;
}
case(G1V2):
{
ret = nc_def_dim(ncid, "longitude", nx, &xdimid);
if(ret != NC_NOERR) return "Can't create x-dimension in the netcdf file";
ret = nc_def_dim(ncid, "latitude", ny, &ydimid);
if(ret != NC_NOERR) return "Can't create y-dimension in the netcdf file";
ret = nc_def_var(ncid, "longitude", NC_FLOAT, 1, &xdimid, &xid);
if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file";
ret = nc_def_var(ncid, "latitude", NC_FLOAT, 1, &ydimid, &yid);
if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file";
break;
}
case(G2V2):
{
ret = nc_def_dim(ncid, "longitude", nx, &xdimid);
if(ret != NC_NOERR) return "Can't create x-dimension in the netcdf file";
ret = nc_def_dim(ncid, "latitude", ny, &ydimid);
if(ret != NC_NOERR) return "Can't create y-dimension in the netcdf file";
ret = nc_def_var(ncid, "longitude", NC_FLOAT, 2, dimid, &xid);
if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file";
ret = nc_def_var(ncid, "latitude", NC_FLOAT, 2, dimid, &yid);
if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file";
break;
}
case(UNKNOWN): return "Can't determine file type";
}
ret = nc_def_var_deflate(ncid, xid, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for longitude variable in the netcdf file";
ret = nc_def_var_deflate(ncid, yid, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for latitude variable in the netcdf file";
text = "longitude";
ret = nc_put_att_text(ncid, xid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of longitude variable in the netcdf file";
text = "latitude";
ret = nc_put_att_text(ncid, yid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of latitude variable in the netcdf file";
type = stype;
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(HaveVar(name)) return "Variable " + name + " already defined";
struct Var v(name, 0);
int ret;
if(type == G1V1)
ret = nc_def_var(ncid, v.name.Buf(), NC_FLOAT, 1, &xdimid, &v.id);
else
ret = nc_def_var(ncid, v.name.Buf(), NC_FLOAT, 2, dimid, &v.id);
if(ret != NC_NOERR) return "Can't create " + v.name + " variable in the netcdf file";
ret = nc_def_var_deflate(ncid, v.id, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for " + v.name + " variable in the netcdf file";
if(stname.Exist()) ret = nc_put_att_text(ncid, v.id, "standard_name", stname.Len() + 1, stname.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of " + v.name + " variable in the netcdf file";
if(lname.Exist()) ret = nc_put_att_text(ncid, v.id, "long_name", lname.Len() + 1, lname.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of " + v.name + " variable in the netcdf file";
if(units.Exist()) ret = nc_put_att_text(ncid, v.id, "units", units.Len() + 1, units.Buf());
if(ret != NC_NOERR) return "Can't write units attribute of " + v.name + " variable in the netcdf file";
if(comment.Exist()) ret = nc_put_att_text(ncid, v.id, "comment", comment.Len() + 1, comment.Buf());
if(ret != NC_NOERR) return "Can't write comment attribute of " + v.name + " variable in the netcdf file";
ret = nc_put_att_float(ncid, v.id, "_FillValue", NC_FLOAT, 1, &fill);
if(ret != NC_NOERR) return "Can't write _FillValue attribute of " + v.name + " variable in the netcdf file";
vars.push_back(v);
return "";
}
Loading…
Cancel
Save