diff --git a/actions/actionint.h b/actions/actionint.h new file mode 100644 index 0000000..ccca732 --- /dev/null +++ b/actions/actionint.h @@ -0,0 +1,210 @@ +#pragma once +#include "BFileR.h" +#include "BFileW.h" +#include "actiondep.h" +#include "interpolation.h" +#include "ncfilew.h" +#include "ncfuncs.h" +#include + +template +concept ReadIsInterpolable = requires { + { + std::declval>().GridPos(0.0, 0.0) + } -> std::convertible_to; +}; + +ADD_ACTION(INT, int, (ReadPSupported || ReadSupported)&&ReadIsInterpolable); + +template MString ActionINT::DoAction(const CLArgs& args, D& ds) +{ + struct TPoint + { + struct Point2D p; + MDateTime t; + real lon0, lat0, dtime; + }; + + MDateTime refdate("1980-01-01"); + + if(!args.contains("input")) return "No input file given"; + michlib::BFileR in; + if(in.Open(args.at("input")) != ERR_NOERR) return "Can't open input file " + args.at("input"); + std::vector points; + struct Region reg; + + { + michlib::uint loncol = in.Columns(), latcol = in.Columns(), timecol = in.Columns(); + michlib::uint lon0col = in.Columns(), lat0col = in.Columns(), dtimecol = in.Columns(); + for(uint i = 0; i < in.Columns(); i++) + { + if(in.ColumnName(i + 1) == "lon") loncol = i; + if(in.ColumnName(i + 1) == "lat") latcol = i; + if(in.ColumnName(i + 1) == "rtime") timecol = i; + if(in.ColumnName(i + 1) == "lon0") lon0col = i; + if(in.ColumnName(i + 1) == "lat0") lat0col = i; + if(in.ColumnName(i + 1) == "time") dtimecol = i; + } + if(loncol >= in.Columns()) return "Lon column not found in input file " + args.at("input"); + if(latcol >= in.Columns()) return "Lat column not found in input file " + args.at("input"); + if(timecol >= in.Columns()) return "Time column not found in input file " + args.at("input"); + if(lon0col >= in.Columns()) return "Lon0 column not found in input file " + args.at("input"); + if(lat0col >= in.Columns()) return "Lat0 column not found in input file " + args.at("input"); + if(dtimecol >= in.Columns()) return "Age column not found in input file " + args.at("input"); + + points.resize(in.Rows()); + reg.lonb = reg.lone = in[loncol][0]; + reg.latb = reg.late = in[latcol][0]; + for(uint i = 0; i < in.Rows(); i++) + { + points[i].p.x = in[loncol][i]; + points[i].p.y = in[latcol][i]; + points[i].t = refdate + static_cast(michlib::Round(in[timecol][i] * 86400)); + points[i].lon0 = in[lon0col][i]; + points[i].lat0 = in[lat0col][i]; + points[i].dtime = in[dtimecol][i]; + + reg.lonb = std::min(reg.lonb, points[i].p.x); + reg.lone = std::max(reg.lone, points[i].p.x); + reg.latb = std::min(reg.latb, points[i].p.y); + reg.late = std::max(reg.late, points[i].p.y); + } + std::ranges::sort(points, {}, &TPoint::t); + } + in.Close(); + + auto resop = ds.Open(args); + if(resop.Exist()) return "Can't open source: " + resop; + + michlib_internal::ParameterListEx pars; + pars.UsePrefix(""); + pars.SetParameter("source", args.at("source")); + pars.SetParameter("history", args.at("_cmdline")); + + MString varstring; + if(args.contains("var")) + varstring = args.at("var"); + else + { + if(args.contains("vars")) + varstring = args.at("vars"); + else + { + if constexpr(HasDefVars) + varstring = ds.DefaultVars(); + else + return "Variables not specified"; + } + } + + auto vlist = michlib::Split_on_words(varstring, " ,", false); + std::vector vnames{std::make_move_iterator(std::begin(vlist)), std::make_move_iterator(std::end(vlist))}; + { + std::set used; + for(const auto& vname: vnames) + { + if(used.contains(vname)) return "Duplicate variable " + vname + " in list " + varstring; + if(ds.CheckVar(vname) == VarPresence::NONE) return "Variable " + vname + " not exists in this dataset"; + used.insert(vname); + } + } + + pars.SetParameter("variables", varstring); + + //int compress = 3; + //if(args.contains("compress")) compress = args.at("compress").ToInt(); + + std::unique_ptr sourcepars; + if constexpr(ParametersSupported) + { + if constexpr(ParametersRequiredRegion) + { + auto [p, err] = ds.Parameters(pars, args, reg); + if(err.Exist()) return err; + sourcepars.reset(p); + } + else + { + auto [p, err] = ds.Parameters(pars, args); + if(err.Exist()) return err; + sourcepars.reset(p); + } + } + auto p = sourcepars.get(); + + MString name = args.contains("out") ? args.at("out") : "out.bin"; + //MString outfmt = args.contains("format") ? args.at("format") : (GetExt(name) == "nc" ? "nc" : "bin"); + + size_t loind = ds.NTimes(), hiind = ds.NTimes(), curind = 0; + michlib::BFileW fw; + const MDateTime first = ds.Time(0), last = ds.Time(ds.NTimes() - 1); + + std::vector>> lo, hi; + + { + size_t ic = 0; + fw.Create(name, 6 + vnames.size()); + fw.SetColumnName(++ic, "lon0"); + fw.SetColumnName(++ic, "lat0"); + fw.SetColumnName(++ic, "lon"); + fw.SetColumnName(++ic, "lat"); + fw.SetColumnName(++ic, "time"); + fw.SetColumnName(++ic, "rtime"); + for(size_t i = 0; i < vnames.size(); i++) fw.SetColumnName(++ic, vnames[i]); + fw.SetParameters(pars); + } + + for(size_t i = 0; i < points.size(); i++) + { + if(points[i].t < first || points[i].t > last) continue; + while(!(points[i].t >= ds.Time(curind) && points[i].t <= ds.Time(curind + 1))) curind++; + if(curind != loind && curind != hiind) + { + loind = curind; + hiind = curind + 1; + { + auto temp = Read(ds, vnames, p, loind); + if(temp.size() != vnames.size()) return "Can't read data"; + //lo.resize(temp.size()); + for(size_t j = 0; j < temp.size(); j++) lo.emplace_back(std::move(temp[j])); + } + { + auto temp = Read(ds, vnames, p, hiind); + if(temp.size() != vnames.size()) return "Can't read data"; + //hi.resize(temp.size()); + for(size_t j = 0; j < temp.size(); j++) hi.emplace_back(std::move(temp[j])); + } + } + else if(curind == hiind) + { + loind = curind; + hiind = curind + 1; + lo = std::move(hi); + { + auto temp = Read(ds, vnames, p, hiind); + if(temp.size() != vnames.size()) return "Can't read data"; + //hi.resize(temp.size()); + for(size_t j = 0; j < temp.size(); j++) hi.emplace_back(std::move(temp[j])); + } + } + auto step = ds.Time(hiind) - ds.Time(loind); + auto delta = points[i].t - ds.Time(loind); + real trel = static_cast(delta) / static_cast(step); + + fw.Write(points[i].lon0); + fw.Write(points[i].lat0); + fw.Write(points[i].p.x); + fw.Write(points[i].p.y); + fw.Write(points[i].dtime); + fw.Write(DeltaDates(refdate, points[i].t)); + for(size_t iv = 0; iv < lo.size(); iv++) + { + real vlo = lo[iv](points[i].p); + real vhi = hi[iv](points[i].p); + fw.Write(vlo + (vhi - vlo) * trel); + } + } + fw.Finalize(); + fw.Close(); + return ""; +}; diff --git a/include/basedata.h b/include/basedata.h index 5b8948b..38c3e98 100644 --- a/include/basedata.h +++ b/include/basedata.h @@ -15,6 +15,19 @@ struct Region real lonb, lone, latb, late; }; +struct GridPoint +{ + size_t ix = 0, iy = 0; + real x = -1.0, y = -1.0; + + bool Valid() const { return x >= 0.0 && y >= 0.0; } +}; + +struct Point2D +{ + real x, y; +}; + enum class VarPresence { NONE, diff --git a/include/interpolation.h b/include/interpolation.h new file mode 100644 index 0000000..c544925 --- /dev/null +++ b/include/interpolation.h @@ -0,0 +1,51 @@ +#pragma once +#include "basedata.h" +#include "comdefs.h" + +using michlib::real; + +template class LinearInterpolator: public D +{ + public: + LinearInterpolator(D&& d): D(std::move(d)) {} + + real operator()(const struct Point2D& in) const + { + const auto gp = D::GridPos(in.x, in.y); + if(!gp.Valid()) return NAN; + + if(gp.x == 0.0 && gp.y == 0.0) return D::V(gp.ix, gp.iy); + + real v00; + real v10; + real v01; + real v11; + bool isfill = false; + size_t fx, fy; + + // Count fills + for(size_t ix = 0; ix <= 1; ix++) + for(size_t iy = 0; iy <= 1; iy++) + { + if(isfill && D::IsFill(gp.ix + ix, gp.iy + iy)) return NAN; + if(D::IsFill(gp.ix + ix, gp.iy + iy)) + { + fx = ix; + fy = iy; + isfill = true; + } + } + + v00 = D::V(gp.ix, gp.iy); + v10 = D::V(gp.ix + 1, gp.iy); + v01 = D::V(gp.ix, gp.iy + 1); + v11 = D::V(gp.ix + 1, gp.iy + 1); + + if(isfill && fx == 0 && fy == 0) v00 = (v10 + v01 + v11) / 3.0; + if(isfill && fx == 1 && fy == 0) v10 = (v00 + v01 + v11) / 3.0; + if(isfill && fx == 0 && fy == 1) v01 = (v10 + v00 + v11) / 3.0; + if(isfill && fx == 1 && fy == 1) v11 = (v10 + v01 + v00) / 3.0; + + return v00 + gp.y * (v01 - v00) + gp.x * ((v10 - v00) + gp.y * (v00 - v01 + v11 - v10)); + }; +}; diff --git a/include/simple2ddata.h b/include/simple2ddata.h index 9347ec4..b4d33f2 100644 --- a/include/simple2ddata.h +++ b/include/simple2ddata.h @@ -42,6 +42,23 @@ class Simple2DData: public BaseData real XStep() const { return xstep; } real YStep() const { return ystep; } + + struct GridPoint GridPos(real x, real y) const + { + struct GridPoint out; + + if(!*this) return out; + if(x < x0 || x > x0 + (nx - 1) * xstep || y < y0 || y > y0 + (ny - 1) * ystep) return out; + + real rx = x - x0; + real ry = y - y0; + out.ix = static_cast(michlib::Floor(rx / xstep)); + out.iy = static_cast(michlib::Floor(ry / ystep)); + out.x = rx - out.ix * xstep; + out.y = ry - out.iy * ystep; + + return out; + } }; class Rect2DData: public BaseData