Begin work on interpolation code

actions/actioninterpolate2d.cpp

@@ -0,0 +1,10 @@
+#define MICHLIB_NOSOURCE
+#include "actioninterpolate2d.h"
+
+InterpolateMethods::Mode InterpolateMethods::GetMode(const CLArgs& args)
+{
+ MString mode = args.contains("mode") ? args.at("mode") : "vylet_start";
+ if(mode == "vylet_start") return Mode::VYLETSTART;
+ if(mode == "vylet_end") return Mode::VYLETEND;
+ return Mode::INVALID;
+}

actions/actioninterpolate2d.h

@@ -0,0 +1,181 @@
+#pragma once
+#include "BFileR.h"
+#include "actiondep.h"
+
+//using michlib::message;
+
+class InterpolateMethods
+{
+ protected:
+ enum class Mode
+ {
+  INVALID,
+  VYLETSTART,
+  VYLETEND
+ };
+ struct DataPoint
+ {
+  real      lon, lat;
+  MDateTime t;
+ };
+
+ class VyletFile
+ {
+  michlib::BFileR         fr;
+  bool                    invtime;
+  bool                    usestart;
+  MDateTime               beg, end, t0;
+  michlib::CompiledParser xy2lon, xy2lat;
+  real                    x, y;
+
+  public:
+  MString Open(const MString& name, Mode mode);
+ };
+
+ static Mode GetMode(const CLArgs& args);
+};
+
+ADD_ACTION(Interpolate, interpolate, CanInterpolate2D<Source>, InterpolateMethods);
+
+template<class D> MString ActionInterpolate::DoAction(const CLArgs& args, D& ds)
+{
+ Mode mode = GetMode(args);
+ if(mode == Mode::INVALID) return "Unknown mode";
+
+ if(!args.contains("in")) return "No vylet data file specified";
+ MString in = args.at("in");
+
+ if(!args.contains("out")) return "No output file specified";
+ MString out = args.at("out");
+
+ michlib::BFileR fr;
+ if(fr.Open(in) != ERR_NOERR) return "Can't open file " + in;
+
+ std::vector<struct DataPoint> datapoints;
+ {
+  bool invtime;
+  {
+   fr.UsePrefix("");
+   MString method = fr.ParameterSValue("Method", "");
+   if(method == "Bicubic" || method == "BicubicL")
+    invtime = false;
+   else if(method == "BicubicI" || method == "BicubicIL")
+    invtime = true;
+   else
+    return "Unknown method in the file " + in;
+  }
+
+  MDateTime beg, end;
+  {
+   fr.UsePrefix("Datafile_Info");
+   MString s;
+   s = fr.ParameterSValue("BeginDate", "");
+   beg.FromString(s);
+   s = fr.ParameterSValue("EndDate", "");
+   end.FromString(s);
+  }
+
+  MDateTime t0;
+  {
+   fr.UsePrefix("");
+   auto tbeg = static_cast<time_t>(fr.ParameterRValue("tbeg", 0.0) * 86400);
+   if(invtime)
+   {
+    t0 = end;
+    t0.AddSeconds(-tbeg);
+   }
+   else
+   {
+    t0 = beg;
+    t0.AddSeconds(tbeg);
+   }
+  }
+
+  michlib::CompiledParser xy2lon, xy2lat;
+  real                    x, y;
+  michlib::ParserVars     pv;
+
+  pv["x"] = &x;
+  pv["y"] = &y;
+
+  fr.UsePrefix("Datafile_Info");
+  if(!ArifmeticCompiler(fr.ParameterSValue("xy2lon", ""), xy2lon, &pv)) return "Can't find xy2lon in the file " + in;
+  if(!ArifmeticCompiler(fr.ParameterSValue("xy2lat", ""), xy2lat, &pv)) return "Can't find xy2lat in the file " + in;
+
+  auto resop = ds.Open(args);
+  if(resop.Exist()) return "Can't open source: " + resop;
+
+  for(size_t i = 0; i < fr.Rows(); i++)
+  {
+   if(mode == Mode::VYLETSTART)
+   {
+    x = fr[0][i];
+    y = fr[1][i];
+    real lon, lat;
+    xy2lon.Run(lon);
+    xy2lat.Run(lat);
+    datapoints.emplace_back(lon, lat, t0);
+   }
+   if(mode == Mode::VYLETEND)
+   {
+    x = fr[2][i];
+    y = fr[3][i];
+    real lon, lat;
+    xy2lon.Run(lon);
+    xy2lat.Run(lat);
+    time_t    t    = static_cast<time_t>(fr[4][i] * 86400);
+    MDateTime time = invtime ? end : beg;
+    time.AddSeconds((invtime ? -1 : 1) * t);
+    datapoints.emplace_back(lon, lat, time);
+   }
+  }
+ }
+
+ if(datapoints.size() == 0) return "No data readed from file " + in;
+
+ struct Region reg;
+ reg.lonb = datapoints[0].lon;
+ reg.lone = datapoints[0].lon;
+ reg.latb = datapoints[0].lat;
+ reg.late = datapoints[0].lat;
+
+ for(const auto& p: datapoints)
+ {
+  if(p.lon < reg.lonb) reg.lonb = p.lon;
+  if(p.lat < reg.latb) reg.latb = p.lat;
+  if(p.lon > reg.lone) reg.lone = p.lon;
+  if(p.lat > reg.late) reg.late = p.lat;
+ }
+
+ return "";
+ /*
+
+ 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);
+
+ std::unique_ptr<const BaseParameters> sourcepars;
+ if constexpr(ParametersSupported<D>)
+ {
+  if constexpr(ParametersRequiredRegion<D>)
+  {
+   auto [p, err] = ds.Parameters(pars, args, reg);
+   if(err.Exist()) return err;
+   sourcepars.reset(p);
+  }
+  else
+  {
+   auto [p, err] = ds.Parameters(pars, args);
+   if(err.Exist()) return err;
+   sourcepars.reset(p);
+  }
+ }
+ auto p = sourcepars.get();
+
+ auto data = Read(ds, vname, p, tindexes);
+ if(!data) return "Can't read data";
+
+ return "";
+ */
+};

include/basedata.h

@@ -15,6 +15,12 @@ struct Region
  real lonb, lone, latb, late;
 };
 
+struct GridPointLocation
+{
+ size_t ix, iy;
+ real   x, y;
+};
+
 class BaseData
 {
  protected:

include/simple2ddata.h

@@ -1,6 +1,8 @@
 #pragma once
 #include "basedata.h"
 
+using michlib::Floor;
+
 class Simple2DData: public BaseData
 {
  real   x0 = 0.0, y0 = 0.0;
@@ -42,6 +44,16 @@ class Simple2DData: public BaseData
 
  real XStep() const { return xstep; }
  real YStep() const { return ystep; }
+
+ auto GridPoint(real lon, real lat) const
+ {
+  struct GridPointLocation loc;
+  loc.ix = static_cast<size_t>(Floor((lon - x0) / xstep));
+  loc.iy = static_cast<size_t>(Floor((lat - y0) / ystep));
+  loc.x  = lon - x0 - loc.ix * xstep;
+  loc.y  = lat - y0 - loc.iy * ystep;
+  return loc;
+ }
 };
 
 class Rect2DData: public BaseData

include/traits.h

@@ -167,3 +167,10 @@ concept ReadIs1DGeoArray = requires {
   std::declval<ReadType<T>>().Lat(0)
  } -> std::convertible_to<real>;
 };
+
+template<class T>
+concept CanInterpolate2D = requires {
+ {
+  std::declval<ReadType<T>>().GridPoint(0.0, 0.0)
+ } -> std::convertible_to<struct GridPointLocation>;
+};