odm/include/Data3D.h

#pragma once
#include "comdefs.h"
#include "udunitscpp.h"
#include <ranges>
#include <variant>

using michlib::real;

class ProjectionInternals
{
 public:
 enum class Type
 {
  IRREGULAR,   // Набор точек
  IRREGULARXY, // Прямоугольная сетка с неравномерным шагом
  EQC          // Прямоугольная сетка с равномерным шагом
 };

 enum class Metric
 {
  SPHERIC,
  PLANE
 };

 template<Type type> class ProjData;

 private:
 template<class C> struct PrTypeHold;

 template<Type t> struct PrTypeHold<ProjData<t>>
 {
  static constexpr Type type = t;
 };

 public:
 class Sizes
 {
  size_t nx, ny;

  Sizes() = delete;

  public:
  Sizes(size_t nx_, size_t ny_): nx(nx_), ny(ny_) {}

  size_t Nx() const { return nx; }
  size_t Ny() const { return ny; }
  size_t N() const { return nx * ny; }
 };

 template<Type type> static constexpr Metric metric = Metric::SPHERIC;

 // Флаги
 // orthogonal - ортогональная сетка
 // rectangular - прямоугольная сетка, широта и долгота зависят от одного индекса
 // geostep - шаги сетки по широте и долготе фиксированы
 // analytic - существует непрерывное преобразование сеточных координат в широту и долготу и наоборот

 template<Type type> static constexpr bool orthogonal = true;
 template<Type type>
 static constexpr bool rectangular = requires(ProjData<type> pr, size_t i) {
  { pr.Lon(i) } -> std::same_as<real>;
  { pr.Lat(i) } -> std::same_as<real>;
 };
 template<Type type>
 static constexpr bool geostep = requires(ProjData<type> pr) {
  { pr.LonStep() } -> std::same_as<real>;
  { pr.LatStep() } -> std::same_as<real>;
 };
 template<Type type>
 static constexpr bool analytic = requires(ProjData<type> pr, real x, real y, real lon, real lat) {
  { pr.LonR(x, y) } -> std::same_as<real>;
  { pr.LatR(x, y) } -> std::same_as<real>;
  { pr.X(lon, lat) } -> std::same_as<real>;
  { pr.Y(lon, lat) } -> std::same_as<real>;
 };

 template<Type type, class... Args>
 static constexpr bool validargs = requires(Args... args) {
  { ProjData<type>(args...) };
 };

 template<class C> static constexpr Type ProjType = PrTypeHold<C>::type;

 using Projection = std::variant<ProjData<Type::IRREGULAR>, ProjData<Type::IRREGULARXY>, ProjData<Type::EQC>>;
};

template<> constexpr bool ProjectionInternals::orthogonal<ProjectionInternals::Type::IRREGULAR> = false;

template<> class ProjectionInternals::ProjData<ProjectionInternals::Type::IRREGULAR>: public ProjectionInternals::Sizes
{
 std::vector<real> datalon, datalat;

 ProjData() = delete;

 public:
 ProjData(size_t nx, size_t ny, const std::vector<real>& lons, const std::vector<real>& lats): Sizes{nx, ny}, datalon(lons), datalat(lats) {}
 ProjData(size_t nx, size_t ny, std::vector<real>&& lons, std::vector<real>&& lats): Sizes{nx, ny}, datalon(std::move(lons)), datalat(std::move(lats)) {}

 real Lon(size_t ix, size_t iy) const { return datalon[Nx() * iy + ix]; }
 real Lat(size_t ix, size_t iy) const { return datalat[Nx() * iy + ix]; }
};

template<> class ProjectionInternals::ProjData<ProjectionInternals::Type::IRREGULARXY>: public ProjectionInternals::Sizes
{
 std::vector<real> datalon, datalat;

 ProjData() = delete;

 public:
 ProjData(const std::vector<real>& lons, const std::vector<real>& lats): Sizes{lons.size(), lats.size()}, datalon(lons), datalat(lats) {}
 ProjData(std::vector<real>&& lons, std::vector<real>&& lats): Sizes{lons.size(), lats.size()}, datalon(std::move(lons)), datalat(std::move(lats)) {}

 real Lon(size_t ix) const { return datalon[ix]; }
 real Lat(size_t iy) const { return datalat[iy]; }
};

template<> class ProjectionInternals::ProjData<ProjectionInternals::Type::EQC>: public ProjectionInternals::Sizes
{
 real lon0, lat0, lonstep, latstep;

 ProjData() = delete;

 public:
 ProjData(size_t nx, size_t ny, real lon0_, real lat0_, real lonstep_, real latstep_): Sizes{nx, ny}, lon0(lon0_), lat0(lat0_), lonstep(lonstep_), latstep(latstep_) {}
 ProjData(std::ranges::sized_range auto&& lons, std::ranges::sized_range auto&& lats):
     Sizes{lons.size(), lats.size()},
     lon0(lons[0]),
     lat0(lats[0]),
     lonstep((lons.back() - lons.front()) / (lons.size() - 1)),
     latstep((lats.back() - lats.front()) / (lats.size() - 1))
 {
 }

 real Lon(size_t ix) const { return lon0 + ix * lonstep; }
 real Lat(size_t iy) const { return lat0 + iy * latstep; }

 real LonStep() const { return lonstep; }
 real LatStep() const { return latstep; }

 real LonR(real x, [[maybe_unused]] real y) const { return lon0 + x * lonstep; }
 real LatR([[maybe_unused]] real x, real y) const { return lat0 + y * latstep; }

 real X(real lon, [[maybe_unused]] real lat) const { return (lon - lon0) / lonstep; }
 real Y([[maybe_unused]] real lon, real lat) const { return (lat - lat0) / latstep; }
};

class Projection: public ProjectionInternals, public ProjectionInternals::Projection
{
 Projection()                  = delete;
 Projection(const Projection&) = delete;

 template<Type t> Projection(ProjData<t>&& pr): ProjectionInternals(), ProjectionInternals::Projection(std::move(pr)) {}

 public:
 Projection(Projection&&) = default;

 template<Type t, class... Args> static Projection Create(Args&&... args)
 {
  if constexpr(t == Type::IRREGULAR)
  {
   if constexpr(validargs<Type::IRREGULAR, Args...>)
    return Projection(ProjData<Type::IRREGULAR>(std::forward<Args>(args)...));
   else
    static_assert(false, "Incorrect call of Projection::Create");
  }
  else if constexpr(t == Type::IRREGULARXY)
  {
   if constexpr(validargs<Type::IRREGULARXY, Args...>)
    return Projection(ProjData<Type::IRREGULARXY>(std::forward<Args>(args)...));
   else
    static_assert(false, "Incorrect call of Projection::Create");
  }
  else if constexpr(t == Type::EQC)
  {
   if constexpr(validargs<Type::EQC, Args...>)
    return Projection(ProjData<Type::EQC>(std::forward<Args>(args)...));
   else
    static_assert(false, "Incorrect call of Projection::Create");
  }
  else
   static_assert(false, "Incorrect call of Projection::Create");
 };

 // From Sizes
 size_t Nx() const
 {
  return std::visit([](const auto& arg) { return arg.Nx(); }, *this);
 }
 size_t Ny() const
 {
  return std::visit([](const auto& arg) { return arg.Ny(); }, *this);
 }
 size_t N() const
 {
  return std::visit([](const auto& arg) { return arg.N(); }, *this);
 }

 // Lon and Lat (ix,iy)
 real Lon(size_t ix, size_t iy) const
 {
  return std::visit(
      [ix = ix, iy = iy](const auto& arg)
      {
       using T             = std::decay_t<decltype(arg)>;
       constexpr bool rect = Projection::rectangular<ProjType<T>>;
       if constexpr(rect)
        return arg.Lon(ix);
       else
        return arg.Lon(ix, iy);
      },
      *this);
 }
 real Lat(size_t ix, size_t iy) const
 {
  return std::visit(
      [ix = ix, iy = iy](const auto& arg)
      {
       using T             = std::decay_t<decltype(arg)>;
       constexpr bool rect = Projection::rectangular<ProjType<T>>;
       if constexpr(rect)
        return arg.Lat(iy);
       else
        return arg.Lat(ix, iy);
      },
      *this);
 }

 // Lon and Lat (i)
 real Lon(size_t ix) const
 {
  return std::visit(
      [ix = ix](const auto& arg) -> real
      {
       using T             = std::decay_t<decltype(arg)>;
       constexpr bool rect = Projection::rectangular<ProjType<T>>;
       if constexpr(rect)
        return arg.Lon(ix);
       else
        return NAN;
      },
      *this);
 }
 real Lat(size_t iy) const
 {
  return std::visit(
      [iy = iy](const auto& arg) -> real
      {
       using T             = std::decay_t<decltype(arg)>;
       constexpr bool rect = Projection::rectangular<ProjType<T>>;
       if constexpr(rect)
        return arg.Lat(iy);
       else
        return NAN;
      },
      *this);
 }

 // LonR, LatR
 real LonR(real x, real y) const
 {
  return std::visit(
      [x = x, y = y](const auto& arg) -> real
      {
       using T                 = std::decay_t<decltype(arg)>;
       constexpr bool analytic = Projection::analytic<ProjType<T>>;
       if constexpr(analytic)
        return arg.LonR(x, y);
       else
        return NAN;
      },
      *this);
 }
 real LatR(real x, real y) const
 {
  return std::visit(
      [x = x, y = y](const auto& arg) -> real
      {
       using T                 = std::decay_t<decltype(arg)>;
       constexpr bool analytic = Projection::analytic<ProjType<T>>;
       if constexpr(analytic)
        return arg.LatR(x, y);
       else
        return NAN;
      },
      *this);
 }

 // X, Y
 real X(real lon, real lat) const
 {
  return std::visit(
      [lon = lon, lat = lat](const auto& arg) -> real
      {
       using T                 = std::decay_t<decltype(arg)>;
       constexpr bool analytic = Projection::analytic<ProjType<T>>;
       if constexpr(analytic)
        return arg.X(lon, lat);
       else
        return NAN;
      },
      *this);
 }
 real Y(real lon, real lat) const
 {
  return std::visit(
      [lon = lon, lat = lat](const auto& arg) -> real
      {
       using T                 = std::decay_t<decltype(arg)>;
       constexpr bool analytic = Projection::analytic<ProjType<T>>;
       if constexpr(analytic)
        return arg.Y(lon, lat);
       else
        return NAN;
      },
      *this);
 }

 // Type
 Type Type() const
 {
  return std::visit([](const auto& arg) { return ProjType<std::decay_t<decltype(arg)>>; }, *this);
 }

 // Flags
 bool IsOrthogonal() const
 {
  return std::visit([](const auto& arg) { return Projection::orthogonal<ProjType<std::decay_t<decltype(arg)>>>; }, *this);
 }
 bool IsRectangular() const
 {
  return std::visit([](const auto& arg) { return Projection::rectangular<ProjType<std::decay_t<decltype(arg)>>>; }, *this);
 }
 bool IsAnalytic() const
 {
  return std::visit([](const auto& arg) { return Projection::analytic<ProjType<std::decay_t<decltype(arg)>>>; }, *this);
 }

 // Metric
 auto Metric() const
 {
  return std::visit([](const auto& arg) { return Projection::metric<ProjType<std::decay_t<decltype(arg)>>>; }, *this);
 }
};

class VerticalInternals
{
 public:
 enum class Type
 {
  LAYERS,    // Глубина слоя зависит от координат
  HORIZONTS, // Глубина слоя не зависит от координат
  GRID,      // Равномерная сетка глубин
  FIXED,     // Один уровень фиксированой глубины
  UNDEFINED  // Глубина бессмысленна, например, для ssh или толщины перемешанного слоя
 };

 template<Type type> class VertData;

 private:
 template<class C> struct VTypeHold;

 template<Type t> struct VTypeHold<VertData<t>>
 {
  static constexpr Type type = t;
 };

 public:
 class Sizes
 {
  size_t nz;

  Sizes() = delete;

  public:
  Sizes(size_t nz_): nz(nz_) {}

  size_t Nz() const { return nz; }
 };

 // Флаги
 // up2down - слои заглубляются с увеличением номера
 // isotropic - глубина слоя не зависит от горизонтальных координат
 // zstep - постоянный шаг между слоями
 // valid - глубина определена

 template<Type type> static constexpr bool up2down = true;
 template<Type type>
 static constexpr bool isotropic = requires(VertData<type> prz, size_t i) {
  { prz.Depth(i) } -> std::same_as<real>;
 };
 template<Type type>
 static constexpr bool zstep = requires(VertData<type> prz) {
  { prz.ZStep() } -> std::same_as<real>;
 };
 template<Type type>
 static constexpr bool valid = requires(VertData<type> prz, size_t ix, size_t iy, size_t iz) {
  { prz.Depth(ix, iy, iz) } -> std::same_as<real>;
 } || isotropic<type>;

 template<Type type, class... Args>
 static constexpr bool validargs = requires(Args... args) {
  { VertData<type>(args...) };
 };

 template<class C> static constexpr Type VertType = VTypeHold<C>::type;

 using Vertical = std::variant<VertData<Type::LAYERS>, VertData<Type::HORIZONTS>, VertData<Type::GRID>, VertData<Type::FIXED>, VertData<Type::UNDEFINED>>;
};

template<> class VerticalInternals::VertData<VerticalInternals::Type::LAYERS>: public VerticalInternals::Sizes
{
 size_t            nx, ny;
 std::vector<real> zdata;

 VertData() = delete;

 public:
 VertData(size_t nx_, size_t ny_, size_t nz_, const std::vector<real>& zdata_): Sizes(nz_), nx(nx_), ny(ny_), zdata(zdata_) {}
 VertData(size_t nx_, size_t ny_, size_t nz_, std::vector<real>&& zdata_): Sizes(nz_), nx(nx_), ny(ny_), zdata(std::move(zdata_)) {}

 real Depth(size_t ix, size_t iy, size_t iz) const { return zdata[nx * ny * iz + nx * iy + ix]; }
};

template<> class VerticalInternals::VertData<VerticalInternals::Type::HORIZONTS>: public VerticalInternals::Sizes
{
 std::vector<real> zdata;

 VertData() = delete;

 public:
 VertData(const std::vector<real>& zdata_): Sizes(zdata_.size()), zdata(zdata_) {}
 VertData(std::vector<real>&& zdata_): Sizes(zdata_.size()), zdata(std::move(zdata_)) {}

 real Depth(size_t iz) const { return zdata[iz]; }
};

template<> class VerticalInternals::VertData<VerticalInternals::Type::GRID>: public VerticalInternals::Sizes
{
 real z0, zstep;

 VertData() = delete;

 public:
 VertData(real z0_, real zstep_, size_t nz_): Sizes(nz_), z0(z0_), zstep(zstep_) {}
 VertData(real zstep_, size_t nz_): VertData(0.0, zstep_, nz_) {}

 real Depth(size_t iz) const { return z0 + iz * zstep; }
 real ZStep() const { return zstep; }
};

template<> class VerticalInternals::VertData<VerticalInternals::Type::FIXED>
{
 real z;

 VertData() = delete;

 public:
 VertData(real z_): z(z_) {}

 real Depth([[maybe_unused]] size_t iz) const { return z; }

 constexpr size_t Nz() const { return 1; }
};

template<> class VerticalInternals::VertData<VerticalInternals::Type::UNDEFINED>
{
 public:
 VertData() = default;

 //real Depth([[maybe_unused]] size_t iz) const { return NAN; }

 constexpr size_t Nz() const { return 1; }
};

class Vertical: public VerticalInternals, public VerticalInternals::Vertical
{
 Vertical()                = delete;
 Vertical(const Vertical&) = delete;

 template<Type t> Vertical(VertData<t>&& prz): VerticalInternals(), VerticalInternals::Vertical(std::move(prz)) {}

 public:
 template<Type t, class... Args> static Vertical Create(Args&&... args)
 {
  if constexpr(t == Type::LAYERS)
  {
   if constexpr(validargs<Type::LAYERS, Args...>)
    return Vertical(VertData<Type::LAYERS>(std::forward<Args>(args)...));
   else
    static_assert(false, "Incorrect call of Vertical::Create");
  }
  else if constexpr(t == Type::HORIZONTS)
  {
   if constexpr(validargs<Type::HORIZONTS, Args...>)
    return Vertical(VertData<Type::HORIZONTS>(std::forward<Args>(args)...));
   else
    static_assert(false, "Incorrect call of Vertical::Create");
  }
  else if constexpr(t == Type::GRID)
  {
   if constexpr(validargs<Type::GRID, Args...>)
    return Vertical(VertData<Type::GRID>(std::forward<Args>(args)...));
   else
    static_assert(false, "Incorrect call of Vertical::Create");
  }
  else if constexpr(t == Type::FIXED)
  {
   if constexpr(validargs<Type::FIXED, Args...>)
    return Vertical(VertData<Type::FIXED>(std::forward<Args>(args)...));
   else
    static_assert(false, "Incorrect call of Vertical::Create");
  }
  else if constexpr(t == Type::UNDEFINED)
  {
   if constexpr(validargs<Type::UNDEFINED, Args...>)
    return Vertical(VertData<Type::UNDEFINED>(std::forward<Args>(args)...));
   else
    static_assert(false, "Incorrect call of Vertical::Create");
  }
  else
   static_assert(false, "Incorrect call of Vertical::Create");
 };

 // From Sizes
 size_t Nz() const
 {
  return std::visit([](const auto& arg) { return arg.Nz(); }, *this);
 }

 // Depth (ix,iy,iz)
 real Depth(size_t ix, size_t iy, size_t iz) const
 {
  return std::visit(
      [ix = ix, iy = iy, iz = iz](const auto& arg) -> real
      {
       using T                  = std::decay_t<decltype(arg)>;
       constexpr bool isotropic = Vertical::isotropic<VertType<T>>;
       constexpr bool valid     = Vertical::valid<VertType<T>>;
       if constexpr(!valid)
        return NAN;
       else if constexpr(isotropic)
        return arg.Depth(iz);
       else
        return arg.Depth(ix, iy, iz);
      },
      *this);
 }

 // Depth (iz)
 real Depth(size_t iz) const
 {
  return std::visit(
      [iz = iz](const auto& arg) -> real
      {
       using T                  = std::decay_t<decltype(arg)>;
       constexpr bool isotropic = Vertical::isotropic<VertType<T>>;
       constexpr bool valid     = Vertical::valid<VertType<T>>;
       if constexpr(isotropic && valid)
        return arg.Depth(iz);
       else
        return NAN;
      },
      *this);
 }

 // Type
 Type Type() const
 {
  return std::visit([](const auto& arg) { return VertType<std::decay_t<decltype(arg)>>; }, *this);
 }

 // Flags
 bool IsIsotropic() const
 {
  return std::visit([](const auto& arg) { return Vertical::isotropic<Vertical::VertType<std::decay_t<decltype(arg)>>>; }, *this);
 }

 bool IsValid() const
 {
  return std::visit([](const auto& arg) { return Vertical::valid<Vertical::VertType<std::decay_t<decltype(arg)>>>; }, *this);
 }
};

class DataMeta
{
 MString unit, stname, lname, comment;

 public:
 struct Store
 {
  MString unit = "", stname = "", lname = "", comment = "";
  real    fillval = NAN;
 };

 private:
 struct Store store;

 public:
 DataMeta(struct Store&& meta): store(std::move(meta)) {}

 const MString& Unit() const { return store.unit; }
 const MString& StandardName() const { return store.stname; }
 const MString& LongName() const { return store.lname; }
 const MString& Comment() const { return store.comment; }

 real Fillval() const { return store.fillval; }

 void SetFillval(real f) { store.fillval = f; }
};

class Data3D: public DataMeta
{
 std::vector<real>           data;
 std::shared_ptr<Projection> proj;
 std::shared_ptr<Vertical>   vert;

 Data3D()              = delete;
 Data3D(const Data3D&) = delete;

 size_t I(size_t ix, size_t iy, size_t iz) const { return ix + iy * Nx() + iz * Nx() * Ny(); }

 Data3D(std::vector<real>&& data_, std::shared_ptr<Projection> proj_, std::shared_ptr<Vertical> vert_, DataMeta::Store&& meta_):
     DataMeta(std::move(meta_)), data(std::move(data_)), proj(proj_), vert(vert_)
 {
 }

 public:
 static Data3D Create(std::vector<real>&& data_, std::shared_ptr<Projection> proj_, std::shared_ptr<Vertical> vert_, DataMeta::Store&& meta_)
 {
  return Data3D(std::move(data_), proj_, vert_, std::move(meta_));
 }
 Data3D(std::shared_ptr<Projection> proj_, std::shared_ptr<Vertical> vert_, DataMeta::Store&& meta_):
     DataMeta(std::move(meta_)), data(proj_->Nx() * proj_->Ny() * vert_->Nz()), proj(proj_), vert(vert_)
 {
 }

 const auto& Data() const { return data; }

 bool IsFill(size_t ix, size_t iy, size_t iz) const { return IsFill(I(ix, iy, iz)); }
 bool IsFill(size_t i) const { return data[i] == Fillval(); }

 size_t Nx() const { return proj->Nx(); }
 size_t Ny() const { return proj->Ny(); }
 size_t Nz() const { return vert->Nz(); }

 const real& V(size_t ix, size_t iy, size_t iz) const { return V(I(ix, iy, iz)); }
 real&       V(size_t ix, size_t iy, size_t iz) { return V(I(ix, iy, iz)); }

 const real& V(size_t i) const { return data[i]; }
 real&       V(size_t i) { return data[i]; }

 const real& operator()(size_t ix, size_t iy, size_t iz) const { return V(ix, iy, iz); }
 real&       operator()(size_t ix, size_t iy, size_t iz) { return V(ix, iy, iz); }
 const real& operator()(size_t i) const { return V(i); }
 real&       operator()(size_t i) { return V(i); }

 real Lon(size_t ix, size_t iy) const { return proj->Lon(ix, iy); }
 real Lat(size_t ix, size_t iy) const { return proj->Lat(ix, iy); }
 real Lon(size_t ix) const { return proj->Lon(ix); }
 real Lat(size_t iy) const { return proj->Lat(iy); }

 real LonR(real x, real y) const { return proj->LonR(x, y); }
 real LatR(real x, real y) const { return proj->LatR(x, y); }

 real X(real lon, real lat) const { return proj->X(lon, lat); }
 real Y(real lon, real lat) const { return proj->Y(lon, lat); }

 real Depth(size_t ix, size_t iy, size_t iz) const { return vert->Depth(ix, iy, iz); }
 real Depth(size_t iz) const { return vert->Depth(iz); }

 auto HType() const { return proj->Type(); }
 auto VType() const { return vert->Type(); }

 // Flags
 bool IsOrthogonal() const { return proj->IsOrthogonal(); }
 bool IsRectangular() const { return proj->IsRectangular(); }
 bool IsAnalytic() const { return proj->IsAnalytic(); }
 bool IsIsotropic() const { return vert->IsIsotropic(); }
 bool IsDepthValid() const { return vert->IsValid(); }

 // Metric
 auto Metric() const { return proj->Metric(); }
};

class Data2D: public DataMeta
{
 std::vector<real>           data;
 std::shared_ptr<Projection> proj;

 Data2D()              = delete;
 Data2D(const Data2D&) = delete;

 size_t I(size_t ix, size_t iy) const { return ix + iy * Nx(); }

 Data2D(std::vector<real>&& data_, std::shared_ptr<Projection> proj_, DataMeta::Store&& meta_): DataMeta(std::move(meta_)), data(std::move(data_)), proj(proj_) {}

 public:
 static Data2D Create(std::vector<real>&& data_, std::shared_ptr<Projection> proj_, DataMeta::Store&& meta_) { return Data2D(std::move(data_), proj_, std::move(meta_)); }
 Data2D(std::shared_ptr<Projection> proj_, DataMeta::Store&& meta_): DataMeta(std::move(meta_)), data(proj_->Nx() * proj_->Ny()), proj(proj_) {}

 const auto& Data() const { return data; }

 bool IsFill(size_t ix, size_t iy) const { return IsFill(I(ix, iy)); }
 bool IsFill(size_t i) const { return data[i] == Fillval(); }

 size_t Nx() const { return proj->Nx(); }
 size_t Ny() const { return proj->Ny(); }

 const real& V(size_t ix, size_t iy) const { return V(I(ix, iy)); }
 real&       V(size_t ix, size_t iy) { return V(I(ix, iy)); }

 const real& V(size_t i) const { return data[i]; }
 real&       V(size_t i) { return data[i]; }

 const real& operator()(size_t ix, size_t iy) const { return V(ix, iy); }
 real&       operator()(size_t ix, size_t iy) { return V(ix, iy); }
 const real& operator()(size_t i) const { return V(i); }
 real&       operator()(size_t i) { return V(i); }

 real Lon(size_t ix, size_t iy) const { return proj->Lon(ix, iy); }
 real Lat(size_t ix, size_t iy) const { return proj->Lat(ix, iy); }
 real Lon(size_t ix) const { return proj->Lon(ix); }
 real Lat(size_t iy) const { return proj->Lat(iy); }

 real LonR(real x, real y) const { return proj->LonR(x, y); }
 real LatR(real x, real y) const { return proj->LatR(x, y); }

 real X(real lon, real lat) const { return proj->X(lon, lat); }
 real Y(real lon, real lat) const { return proj->Y(lon, lat); }

 auto HType() const { return proj->Type(); }

 // Flags
 bool IsOrthogonal() const { return proj->IsOrthogonal(); }
 bool IsRectangular() const { return proj->IsRectangular(); }
 bool IsAnalytic() const { return proj->IsAnalytic(); }

 // Metric
 auto Metric() const { return proj->Metric(); }
};