HYCOM support (dataset=analysis not work for now)

include/HYCOM.h

@@ -0,0 +1,404 @@
+#pragma once
+#include "DataAdapters/ncfilealt.h"
+#include "ParseArgs.h"
+#include "mdatetime.h"
+#include "simple2ddata.h"
+#include "specfunc.h"
+#include <memory>
+#include <regex.h>
+#include <set>
+#include <stdlib.h>
+#include <sys/types.h>
+
+using michlib::Ceil;
+using michlib::errmessage;
+using michlib::Floor;
+using michlib::GPL;
+using michlib::int2;
+using michlib::MDateTime;
+using michlib::NCFileA;
+using michlib::ToGeoDomainPos;
+
+class HYCOMData
+{
+ enum Type
+ {
+  TYPE_UNKNOWN,
+  TYPE_REANALYSIS,
+  TYPE_HINDCAST,
+  TYPE_FORECAST
+ };
+
+ struct Parameters: public BaseParameters
+ {
+  size_t  xb, yb, xe, ye, layer;
+  MString varname;
+  virtual ~Parameters() override = default;
+ };
+
+ std::vector<NCFileA>   nc;
+ std::vector<MString>   urls;
+ std::vector<real>      depths;
+ std::vector<MDateTime> times;
+ size_t                 nx, ny;
+ real                   lonb, latb, lone, late;
+ real                   lonstep, latstep;
+ Type                   type = TYPE_UNKNOWN;
+
+ public:
+ using Data = Simple2DData;
+
+ private:
+ template<class DataType> Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const Parameters* p) const
+ {
+  real     unitmul = 1.0;
+  DataType fill;
+  real     offset = 0.0, scale = 1.0;
+
+  {
+   auto a_fill   = f.A<DataType>(name, "_FillValue");
+   auto a_offset_d = f.A<double>(name, "add_offset");
+   auto a_scale_d  = f.A<double>(name, "scale_factor");
+   auto a_offset_f = f.A<float>(name, "add_offset");
+   auto a_scale_f  = f.A<float>(name, "scale_factor");
+   if(!a_fill) return Data();
+   fill = a_fill;
+   if(a_offset_d) offset = a_offset_d;
+   if(a_scale_d) scale = a_scale_d;
+   if(a_offset_f) offset = a_offset_f;
+   if(a_scale_f) scale = a_scale_f;
+  }
+
+  auto unit = f.A<MString>(name, "units");
+  if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0;
+
+  Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep);
+
+  if(p->xb < p->xe)
+  {
+   auto var = nodepth ? f.V<DataType>(name, {"lon", p->xb, p->xe - p->xb + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
+                      : f.V<DataType>(name, {"lon", p->xb, p->xe - p->xb + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
+
+   if(!var) return Data();
+   if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data();
+
+   for(size_t ix = 0; ix < var.DimLen(0); ix++)
+    for(size_t iy = 0; iy < var.DimLen(1); iy++)
+    {
+     DataType v   = var(ix, iy);
+     data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+    }
+  }
+  else
+  {
+   auto var1 = nodepth ? f.V<DataType>(name, {"lon", p->xb}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
+                       : f.V<DataType>(name, {"lon", p->xb}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
+   auto var2 = nodepth ? f.V<DataType>(name, {"lon", 0, p->xe + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
+                       : f.V<DataType>(name, {"lon", 0, p->xe + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
+   if(!(var1 && var2)) return Data();
+   if((var1.DimLen(0) + var2.DimLen(0)) != data.Nx() || var1.DimLen(1) != data.Ny() || var2.DimLen(1) != data.Ny()) return Data();
+   for(size_t ix = 0; ix < var1.DimLen(0); ix++)
+    for(size_t iy = 0; iy < var1.DimLen(1); iy++)
+    {
+     DataType v   = var1(ix, iy);
+     data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+    }
+   for(size_t ix = 0; ix < var2.DimLen(0); ix++)
+    for(size_t iy = 0; iy < var2.DimLen(1); iy++)
+    {
+     DataType v                    = var2(ix, iy);
+     data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
+    }
+  }
+  return data;
+ }
+
+ static MString VName2IName(const MString& vname)
+ {
+  if(vname == "temp") return "water_temp";
+  if(vname == "sal") return "salinity";
+  if(vname == "ssh") return "surf_el";
+  return "";
+ }
+
+ public:
+ Data Read(const BaseParameters* ip, size_t i) const
+ {
+  if(!isOk()) return Data();
+  bool nodepth = false;
+
+  auto    p    = dynamic_cast<const struct Parameters*>(ip);
+  MString name = VName2IName(p->varname);
+  for(const auto& f: nc)
+  {
+   auto head = f.Header();
+   for(const auto& v: head.Variables())
+    if(v.Name() == name)
+    {
+     if(v.Dimensions().size() == 3) nodepth = true;
+     if(v.Type().Id() == NC_SHORT) return ReadVarRaw<int2>(f, name, i, nodepth, p);
+     if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(f, name, i, nodepth, p);
+    }
+  }
+  return Data();
+ }
+
+ bool HasVar(const MString vname) const
+ {
+  MString name = VName2IName(vname);
+  if(name == "") return false;
+
+  for(const auto& f: nc)
+  {
+   auto head = f.Header();
+   for(const auto& v: head.Variables())
+    if(v.Name() == name) return true;
+  }
+  return false;
+ }
+
+ public:
+ HYCOMData() = default;
+
+ // TODO: RetVal
+ MString Open(const CLArgs& args)
+ {
+  MString dataset = args.contains("dataset") ? args.at("dataset") : "Forecast";
+
+  nc.clear();
+
+  if(dataset == "Forecast")
+  {
+   urls = {"https://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0/FMRC/GLBy0.08_930_FMRC_best.ncd"};
+   type = TYPE_FORECAST;
+  }
+  else if(dataset == "Hindcast")
+  {
+   urls = {"https://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0"};
+   type = TYPE_HINDCAST;
+  }
+  else if(dataset == "Reanalysis")
+  {
+   urls.clear();
+   for(size_t i = 1994; i <= 2015; i++) urls.push_back("https://tds.hycom.org/thredds/dodsC/GLBv0.08/expt_53.X/data/" + MString(i));
+   type = TYPE_REANALYSIS;
+  }
+  else
+   return "Unknown dataset: " + dataset;
+
+  nc.resize(urls.size());
+  for(size_t i = 0; i < urls.size(); i++)
+  {
+   nc[i].Reset(urls[i]+"#cache&noprefetch");
+   if(!nc[i])
+   {
+    nc.clear();
+    return "Can't connect to database";
+   }
+  }
+
+  auto rdepths = nc[0].V<double>("depth");
+  if(!rdepths)
+  {
+   nc.clear();
+   return "Can't read depths";
+  }
+  depths.resize(rdepths.DimLen(0));
+  for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i);
+
+  auto timeD = nc[0].V<double>("time");
+  auto timeF = nc[0].V<float>("time");
+  if(!(timeD || timeF))
+  {
+   nc.clear();
+   return "Can't read times";
+  }
+
+  MDateTime refdate;
+  {
+   auto units = nc[0].Attribute<MString>("time", "units");
+   if(!units)
+   {
+    nc.clear();
+    return "Can't read refdate";
+   }
+   MString rstr;
+   auto    words = michlib::Split_on_words(units);
+   auto    ci    = words.begin();
+   if(ci != words.end()) ci++;              // skip "hours"
+   if(ci != words.end()) ci++;              // skip "since"
+   if(ci != words.end()) rstr = *ci;        // Day
+   if(ci != words.end()) ci++;
+   if(ci != words.end()) rstr += " " + *ci; // Hours
+   if(!refdate.FromString(rstr))
+   {
+    nc.clear();
+    return "Can't parse " + rstr + " to refdate";
+   }
+  }
+
+  timeD ? times.resize(timeD.DimLen(0)) : times.resize(timeF.DimLen(0));
+  for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(timeD ? timeD(i) : timeF(i)) * 3600;
+
+  auto nlon = nc[0].D("lon");
+  auto nlat = nc[0].D("lat");
+  if(!(nlon && nlat))
+  {
+   nc.clear();
+   return "Can't get number of longitudes/latitudes";
+  }
+  nx = nlon;
+  ny = nlat;
+
+  auto lons = nc[0].V<double>("lon");
+  auto lats = nc[0].V<double>("lat");
+  if(!(lons && lats))
+  {
+   nc.clear();
+   return "Can't get longitudes/latitudes";
+  }
+  lonb    = lons(0);
+  latb    = lats(0);
+  lone    = lons(nx - 1);
+  late    = lats(ny - 1);
+  lonstep = (lone - lonb) / (nx - 1);
+  latstep = (late - latb) / (ny - 1);
+
+  return "";
+ }
+
+ MString Info() const
+ {
+  if(!isOk()) return "";
+  MString d;
+  for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")";
+
+  std::set<MString> vars;
+  for(const auto& f: nc)
+  {
+   auto head = f.Header();
+   for(const auto& v: head.Variables())
+   {
+    if(v.Name() == "water_temp") vars.emplace("temp");
+    if(v.Name() == "salinity") vars.emplace("sal");
+    if(v.Name() == "surf_el") vars.emplace("ssh");
+   }
+  }
+  MString svars;
+  {
+   bool first = true;
+   for(const auto& v: vars)
+   {
+    svars += (first ? "" : ", ") + v;
+    first = false;
+   }
+  }
+
+  // clang-format off
+  return
+   "Dataset: " + Title() + "\n" +
+   " Begin date: " + Time(0).ToString() + "\n" +
+   " End date: " + Time(NTimes()-1).ToString() + "\n" +
+   " Time step: " + Timestep() + " seconds\n" +
+   " Time moments: " + NTimes() + "\n" +
+   " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" +
+   " Grid: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" +
+   " Depths:" + d + "\n" +
+   " Supported variables: " + svars;
+  // clang-format on
+ }
+
+ std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const
+ {
+  std::unique_ptr<struct Parameters> ppar{new struct Parameters};
+
+  if(!args.contains("var")) return {nullptr, "Variable not specified"};
+  ppar->varname = args.at("var");
+  if(!HasVar(ppar->varname)) return {nullptr, "Variable " + ppar->varname + " not exists in this dataset"};
+  if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger<size_t>();
+  if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"};
+  real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer);
+
+  if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"};
+
+  {
+   real lon1 = ToGeoDomainPos(args.at("lonb").ToReal());
+   real lon2 = ToGeoDomainPos(args.at("lone").ToReal());
+   real lat1 = ToGeoDomainPos(args.at("latb").ToReal());
+   real lat2 = ToGeoDomainPos(args.at("late").ToReal());
+
+   ppar->yb = static_cast<size_t>(Floor((lat1 - latb) / latstep));
+   ppar->ye = static_cast<size_t>(Ceil((lat2 - latb) / latstep));
+   if(ppar->ye > ny - 1) ppar->ye = ny - 1;
+   if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"};
+
+   ppar->xb = static_cast<size_t>(Floor((lon1 - lonb) / lonstep));
+   ppar->xe = static_cast<size_t>(Ceil((lon2 - lonb) / lonstep));
+
+   if(ppar->xb == ppar->xe) return {nullptr, "Lonb must be not equal late"};
+
+   if(depth < 0.0 || depth > depths.back())
+    ppar->layer = (depth < 0.0) ? 0 : (depths.size() - 1);
+   else
+    for(size_t i = 0; i < depths.size() - 1; i++)
+    {
+     if(depth >= depths[i] && depth <= depths[i + 1])
+     {
+      ppar->layer = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1);
+      break;
+     }
+    }
+  }
+
+  pars.SetParameter("variable", ppar->varname);
+  pars.SetParameter("depth", Depth(ppar->layer));
+  pars.SetParameter("layer", ppar->layer);
+  pars.SetParameter("dataset", Title());
+  pars.SetParameter("lonb", Lon(ppar->xb));
+  pars.SetParameter("latb", Lat(ppar->yb));
+  pars.SetParameter("lone", Lon(ppar->xe));
+  pars.SetParameter("late", Lat(ppar->ye));
+
+  return {ppar.release(), ""};
+ }
+
+ 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
+ }
+
+ 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; }
+ real      Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; }
+ time_t    Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
+ MDateTime Time(size_t i) const
+ {
+  if(!isOk() || i >= times.size()) return MDateTime();
+  return times[i];
+ }
+
+ bool isOk() const { return nc.size() > 0; }
+ explicit operator bool() const { return nc.size() > 0; }
+
+ size_t NDepths() const { return depths.size(); }
+ size_t NTimes() const { return times.size(); }
+
+ MString Title() const
+ {
+  switch(type)
+  {
+  case(TYPE_REANALYSIS): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Reanalysis";
+  case(TYPE_HINDCAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Hindcast";
+  case(TYPE_FORECAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Forecast";
+  default: return "No title";
+  }
+ }
+
+ static real Fillval() { return Data::Fillval(); }
+};

include/odm.h

@@ -1,6 +1,7 @@
 #pragma once
 #include "BFileW.h"
 #include "NEMO.h"
+#include "HYCOM.h"
 #include <utility>
 #include <variant>
 
@@ -9,7 +10,7 @@ using michlib::errmessage;
 using michlib::GPL;
 using michlib::message;
 
-using DataVariants = std::variant<NEMOData>;
+using DataVariants = std::variant<NEMOData,HYCOMData>;
 
 template<class T>
 concept InfoSupported = requires(T t) {
@@ -143,6 +144,13 @@ class Data: public DataVariants
    if(res.Exist()) return "Can't open source " + src + ":\n" + res;
    *this = Data(std::move(data));
   }
+  else if(src == "HYCOM")
+  {
+   HYCOMData data;
+   auto     res = data.Open(args);
+   if(res.Exist()) return "Can't open source " + src + ":\n" + res;
+   *this = Data(std::move(data));
+  }
   else
    return "Unknown source: " + src;
   return "";