odm/include/layereddata.h

#pragma once
#include "gsw.h"
#include "ncfuncs.h"
#include "simple2ddata.h"
#include <memory>

using michlib::Ceil;
using michlib::DetGeoDomain;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;

class LayeredData: public NCFuncs
{
 public:
 using Data = Simple2DData;

 private:
 class NC
 {
  MString                url;
  NCFileA                nc;
  std::vector<MDateTime> times;

  public:
  NC(MString&& newurl): url(std::move(newurl)) { nc.Reset(url + "#cache&noprefetch"); }

  MString ReadTimes(const MString& tname)
  {
   if(!nc) return "File not open";
   auto time = nc.VR(tname);
   if(!time) return "Can't read times";

   MDateTime refdate;
   time_t    step = 0;
   {
    auto units = nc.Attribute<MString>(tname, "units");
    if(!units) return "Can't read refdate";
    auto [rd, st, suc] = Refdate(units);
    if(!suc) return "Can't parse " + units.Get() + " to refdate";
    if(st == 0) return "Can't get timestep from string " + units.Get();
    refdate = rd;
    step    = st;
   }

   times.resize(time.DimLen(0));
   for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(time(i)) * step;
   return "";
  }

  const NCFileA* operator->() const { return &nc; }
  explicit operator bool() const { return nc; }
  MDateTime Begin() const { return times.front(); }
  MDateTime End() const { return times.back(); }

  const NCFileA& Get() const { return nc; }

  const std::vector<MDateTime>& Times() const { return times; }
  size_t                        Index(MDateTime tm) const
  {
   if(tm < Begin() || tm > End()) return 0;
   size_t b = 0, e = times.size() - 1;
   if(tm == times[b]) return b + 1;
   if(tm == times[e]) return e + 1;
   while(e - b > 1)
   {
    size_t c = (e + b) / 2;
    if(tm == times[c]) return c + 1;
    if(tm > times[c])
     b = c;
    else
     e = c;
   }
   return 0;
  }
 };
 std::vector<NC>        nc;
 std::vector<real>      depths;
 std::vector<MDateTime> times;
 struct CoordNames      dname;
 real                   lonb, latb, lone, late;
 real                   lonstep, latstep;
 MString                title;

 class EnvVar
 {
  MString name, oldvalue;
  bool    activated, saved;

  public:
  EnvVar(): activated(false) {}
  ~EnvVar() { Deactivate(); }

  void Activate(const MString& var, const MString& val)
  {
   if(activated) Deactivate();
   name         = var;
   char* curval = getenv(name.Buf());
   if(nullptr == curval)
    saved = false;
   else
   {
    oldvalue = curval;
    saved    = true;
   }
   setenv(name.Buf(), val.Buf(), 1);
  }
  void Deactivate()
  {
   if(!activated) return;
   if(saved)
    setenv(name.Buf(), oldvalue.Buf(), 1);
   else
    unsetenv(name.Buf());
   activated = false;
  }
 };

 EnvVar proxy;

 protected:
 struct Parameters: public BaseParameters
 {
  size_t xb, yb, xe, ye, layer;
  virtual ~Parameters() override = default;
 };

 // TODO: RetVal
 MString Open(const MString& dataset);

 void SetTitle(const MString& newtitle) { title = newtitle; }

 public:
 MString Info() const;

 std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const;

 bool Read(const MString& vname, std::map<MString, Data>& cache, const BaseParameters* ip, size_t i) const;

 bool isOk() const { return nc.size() > 0; }

 explicit operator bool() const { return nc.size() > 0; }

 real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; }

 real Depth(const BaseParameters* ip) const { return Depth(dynamic_cast<const struct Parameters*>(ip)->layer); }

 real Lon(size_t ix) const { return isOk() ? (lonb + ix * lonstep) : -1000.0; }

 real Lat(size_t iy) const { return isOk() ? (latb + iy * latstep) : -1000.0; }

 size_t NDepths() const { return depths.size(); }

 size_t NTimes() const { return times.size(); }

 MDateTime Time(size_t i) const
 {
  if(!isOk() || i >= times.size()) return MDateTime();
  return times[i];
 }

 time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }

 MString Title() const { return title; }

 MString Dump(const struct Parameters* ppar) const
 {
  // clang-format off
  return
   "Current settings:\n" + MString() +
   " Longitudes: from " + Lon(ppar->xb) + " (" + ppar->xb + ") to "+ Lon(ppar->xe) + " (" + ppar->xe + ")\n" +
   " Latitudes: from "  + Lat(ppar->yb) + " (" + ppar->yb + ") to "+ Lat(ppar->ye) + " (" + ppar->ye + ")\n" +
   " Depth: layer " + ppar->layer + ", depth " + Depth(ppar->layer) + " m\n";
  // clang-format on
 }

 VarPresence CheckVar(const MString& vname) const
 {
  return NCFuncs::CheckVar(vname, [this](const MString& vn) { return HaveVar(vn); });
 }

 private:
 template<class DataType> Data ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct Parameters* p) const;

 bool HaveVar(const MString& vname) const
 {
  for(size_t i = 0; i < nc.size(); i++)
   if(NCFuncs::HaveVar(nc[i].Get(), vname)) return true;
  return false;
 }

 std::tuple<MString, size_t, size_t> VarNameLoc(const MString vname, MDateTime tm) const
 {
  for(size_t i = 0; i < nc.size(); i++)
  {
   auto head = nc[i]->Header();
   auto tind = nc[i].Index(tm);
   if(tind == 0) continue;
   for(const auto& v: head.Variables())
   {
    auto stname = nc[i]->A<MString>(v.Name(), "standard_name");
    if(!stname) continue;
    if(StName2Name(stname) == vname) return {v.Name(), i, tind - 1};
   }
  }
  return {"", 0, 0};
 }
};