diff --git a/sources/VYLET.cpp b/sources/VYLET.cpp new file mode 100644 index 0000000..2d4e5b6 --- /dev/null +++ b/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& 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(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(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(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(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(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; +} diff --git a/sources/VYLET.h b/sources/VYLET.h new file mode 100644 index 0000000..9c9089a --- /dev/null +++ b/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 +#include + +using michlib::M_PI; +using michlib::MDateTime; +using michlib::real; +using michlib::Round; + +class VYLETData +{ + std::unique_ptr vylet; + MDateTime start, end, ref; + bool invtime; + real maxx, maxy; + std::vector lons, lats; + std::vector elon, elat; + + MDateTime R2Time(real r) const + { + auto sec = static_cast(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 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& 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); +};