odm/include/simple2ddata.h

#pragma once
#include "basedata.h"
#include "geohelpers.h"

using michlib::M_PI;

class Simple2DData: public BaseData
{
 real   x0 = 0.0, y0 = 0.0;
 size_t nx = 0, ny = 0;
 real   xstep = 0.0, ystep = 0.0;

 static real D(real lon1, real lat1, real lon2, real lat2) { return michlib::GCD(lon1 * M_PI / 180.0, lat1 * M_PI / 180.0, lon2 * M_PI / 180.0, lat2 * M_PI / 180.0) * 6371.0; }

 public:
 using BaseData::IsFill;
 using BaseData::V;
 using BaseData::operator();
 Simple2DData() = default;

 Simple2DData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_, MString&& unit = ""):
     BaseData(nx_ * ny_, std::move(unit)), x0(x0_), y0(y0_), nx(nx_), ny(ny_), xstep(xs_), ystep(ys_)
 {
 }

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

 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 nx; }
 size_t Ny() const { return ny; }

 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 x0 + ix * xstep; }
 real Lat([[maybe_unused]] size_t ix, size_t iy) const { return y0 + iy * ystep; }

 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 x0 + ix * xstep; }
 real Iy2Lat(size_t iy) const { return y0 + iy * ystep; }

 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<decltype(out.ix)>(michlib::Floor(rx / xstep));
  out.iy  = static_cast<decltype(out.ix)>(michlib::Floor(ry / ystep));
  out.x   = rx - out.ix * xstep;
  out.y   = ry - out.iy * ystep;

  return out;
 }

 real dVdX(size_t ix, size_t iy) const
 {
  if(IsFill(ix, iy)) return Fillval();
  if(ix == 0 || IsFill(ix - 1, iy))
  {
   if(IsFill(ix + 1, iy))
    return Fillval();
   else
    return (V(ix + 1, iy) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy));
  }
  if(ix == nx - 1 || IsFill(ix + 1, iy))
  {
   if(IsFill(ix - 1, iy))
    return Fillval();
   else
    return (V(ix, iy) - V(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy));
  }
  return 0.5 * ((V(ix + 1, iy) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
                (V(ix, iy) - V(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
 }
 real dVdY(size_t ix, size_t iy) const
 {
  if(IsFill(ix, iy)) return Fillval();

  if(iy == 0 || IsFill(ix, iy - 1))
  {
   if(IsFill(ix, iy + 1))
    return Fillval();
   else
    return (V(ix, iy + 1) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1));
  }
  if(iy == ny - 1 || IsFill(ix, iy + 1))
  {
   if(IsFill(ix, iy - 1))
    return Fillval();
   else
    return (V(ix, iy) - V(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy));
  }
  return 0.5 * ((V(ix, iy + 1) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
                (V(ix, iy) - V(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
 }

 real Grad(size_t ix, size_t iy) const
 {
  real grx = dVdX(ix, iy);
  real gry = dVdY(ix, iy);
  if(grx == Fillval() || gry == Fillval()) return Fillval();
  return michlib::Hypot(grx, gry);
 }

 real Grad(size_t i) const { return Grad(i % Nx(), i / Nx()); }
};

class Rect2DData: public BaseData
{
 std::vector<real> lon, lat;

 static real D(real lon1, real lat1, real lon2, real lat2) { return michlib::GCD(lon1 * M_PI / 180.0, lat1 * M_PI / 180.0, lon2 * M_PI / 180.0, lat2 * M_PI / 180.0) * 6371.0; }

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

 real dVdX(size_t ix, size_t iy) const
 {
  if(IsFill(ix, iy)) return Fillval();
  if(ix == 0 || IsFill(ix - 1, iy))
  {
   if(IsFill(ix + 1, iy))
    return Fillval();
   else
    return (V(ix + 1, iy) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy));
  }
  if(ix == Nx() - 1 || IsFill(ix + 1, iy))
  {
   if(IsFill(ix - 1, iy))
    return Fillval();
   else
    return (V(ix, iy) - V(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy));
  }
  return 0.5 * ((V(ix + 1, iy) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
                (V(ix, iy) - V(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
 }
 real dVdY(size_t ix, size_t iy) const
 {
  if(IsFill(ix, iy)) return Fillval();

  if(iy == 0 || IsFill(ix, iy - 1))
  {
   if(IsFill(ix, iy + 1))
    return Fillval();
   else
    return (V(ix, iy + 1) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1));
  }
  if(iy == Ny() - 1 || IsFill(ix, iy + 1))
  {
   if(IsFill(ix, iy - 1))
    return Fillval();
   else
    return (V(ix, iy) - V(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy));
  }
  return 0.5 * ((V(ix, iy + 1) - V(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
                (V(ix, iy) - V(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
 }

 real Grad(size_t ix, size_t iy) const
 {
  real grx = dVdX(ix, iy);
  real gry = dVdY(ix, iy);
  if(grx == Fillval() || gry == Fillval()) return Fillval();
  return michlib::Hypot(grx, gry);
 }

 real Grad(size_t i) const { return Grad(i % Nx(), i / Nx()); }
};