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

527 lines
14 KiB

#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;
}