Update action genintfile for work with adapter-based sources

actions/actiongenintfile.h

@@ -1,5 +1,6 @@
 #pragma once
 #include "BFileW.h"
+#include "CF.h"
 #include "actiondep.h"
 #include <memory>
 
@@ -7,9 +8,216 @@ using michlib::BFileW;
 using michlib::Cos;
 using michlib::M_PI;
 
-ADD_ACTION(GenIntFile, genintfile, ReadIsGrid<Source>);
+ADD_ACTION(GenIntFile, genintfile, AdapterSupported<Source> || ReadIsGrid<Source>);
 
-template<class D> MString ActionGenIntFile::DoAction(const CLArgs& args, D& ds)
+template<class D> MString ActionGenIntFile_DoAction(const CLArgs& args, D& data);
+
+template<class D> MString ActionGenIntFile::DoAction(const CLArgs& args, D& data) { return ActionGenIntFile_DoAction(args, data); }
+
+template<class D>
+requires AdapterSupported<D>
+MString ActionGenIntFile_DoAction(const CLArgs& args, D& ds)
+{
+ michlib_internal::ParameterListEx pars;
+ pars.UsePrefix("");
+ pars.SetParameter("source", args.at("source"));
+ pars.SetParameter("history", args.at("_cmdline"));
+
+ Adapter ad;
+ {
+  auto ret = ds.GetAdapter(args, pars);
+  if(!ret) return "Error opening source";
+  ad = std::move(ret.Value());
+ }
+
+ // Modules
+ {
+  auto ret = CF::TEOS10Module(args, pars, ad);
+  if(!ret) return "Error adding TEOS-10 derived variables";
+ }
+ {
+  auto ret = CF::GeostrophicModule(args, pars, ad);
+  if(!ret) return "Error adding geostrophic velocity";
+ }
+
+ if(!args.contains("out")) return "Output file name not specified";
+
+ MString uvar, vvar;
+
+ // Get U and V
+ {
+  MString upar = args.contains("u") ? args.at("u") : "auto";
+  MString vpar = args.contains("v") ? args.at("v") : "auto";
+
+  if(upar == "auto")
+  {
+   if(ad.HasVar("u"))
+    uvar = "u";
+   else if(ad.HasVar("ugs"))
+    uvar = "ugs";
+   else
+    return "Can't find u variable";
+  }
+  else
+   uvar = upar;
+
+  if(vpar == "auto")
+  {
+   MString temp;
+   for(size_t i = 0; i < uvar.Len(); i++) temp += uvar[i] == 'u' ? "v" : uvar.SubStr(i + 1, 1);
+   if(temp != uvar && ad.HasVar(temp))
+    vvar = temp;
+   else
+    return "Can't find v variable";
+  }
+  else
+   vvar = vpar;
+ }
+
+ // Check grids
+ {
+  const auto& proj = ad.Vars().at(uvar).proj;
+  if(proj != ad.Vars().at(vvar).proj) return "U and V variables has different grids";
+  if(proj->Type() != Projection::Type::EQC) return "U and V variables has unsupported grid";
+  if((ad.Is2D(uvar) && ad.Is3D(vvar)) || (ad.Is3D(uvar) && ad.Is2D(vvar))) return "U and V variables has different type";
+  if(ad.Is3D(uvar) && (ad.Vars().at(uvar).vert != ad.Vars().at(vvar).vert)) return "U and V variables has different vertical grids";
+  if(ad.Is3D(uvar) && ad.Vars().at(uvar).vert->Nz() != 1) return "U and V variables has multilayer vertical grid";
+ }
+
+ // Units
+ const UtUnit inu(ad.Vars().at(uvar).info->targetunit), inv(ad.Vars().at(vvar).info->targetunit);
+ const UtUnit cms("cm/s");
+ const auto   cnvu = Converter(inu, cms);
+ const auto   cnvv = Converter(inv, cms);
+
+ if(!cnvu) return "Can't convert units from \"" + ad.Vars().at(uvar).info->targetunit + "\" to \"cm\"";
+ if(!cnvv) return "Can't convert units from \"" + ad.Vars().at(vvar).info->targetunit + "\" to \"cm\"";
+
+ const bool needsconvertu = inu != cms;
+ const bool needsconvertv = inv != cms;
+
+ auto OutU = [needsconvert = needsconvertu, &cnv = cnvu](real r) -> real { return needsconvert ? cnv(r) : r; };
+ auto OutV = [needsconvert = needsconvertv, &cnv = cnvv](real r) -> real { return needsconvert ? cnv(r) : r; };
+
+ const auto& ind   = ad.TimeIndexes();
+ const auto  times = ind | std::views::transform([&times = ad.Times()](size_t i) { return times[i]; });
+ if(!CF::CheckUniform(times)) return "Time is not uniform";
+
+ michlib::message("U is \"", uvar, "\", V is \"", vvar, "\"");
+ pars.SetParameter("U", uvar);
+ pars.SetParameter("V", vvar);
+
+ const real fillval = 1e10;
+
+ MString name = args.at("out");
+ BFileW  fw;
+ fw.Create(name, 2);
+ fw.SetColumnName(1, "u");
+ fw.SetColumnName(2, "v");
+
+ const auto&  proj = ad.Vars().at(uvar).proj;
+ const size_t Nx   = proj->Nx();
+ const size_t Ny   = proj->Ny();
+ {
+  real lonb = proj->Lon(0), latb = proj->Lat(0), lone = proj->Lon(Nx - 1), late = proj->Lat(Ny - 1);
+  fw.UsePrefix("");
+  fw.SetParameter("nx", Nx);
+  fw.SetParameter("ny", Ny);
+  fw.SetParameter("dx", (lone - lonb) / (Nx - 1) * 60.0);
+  fw.SetParameter("dy", (late - latb) / (Ny - 1) * 60.0);
+  fw.SetParameter("nt", times.size());
+  fw.SetParameter("dt", (times[1] - times[0]).D());
+  fw.SetParameter("FillValue", fillval);
+  fw.UsePrefix("Info");
+
+  fw.SetParameter("lonb", lonb);
+  fw.SetParameter("latb", latb);
+  fw.SetParameter("lone", lone);
+  fw.SetParameter("late", late);
+  fw.SetParameter("DataSource", ad.Title());
+  fw.SetParameter("xy2lon", MString(lonb) + "+%x/60");
+  fw.SetParameter("xy2lon_rp", MString("POP 60 DIV ") + lonb + " ADD");
+  fw.SetParameter("lonlat2x", MString("(%lon-") + lonb + ")*60");
+  fw.SetParameter("xy2lat", MString(latb) + "+%y/60");
+  fw.SetParameter("xy2lat_rp", MString("EXCH POP 60 DIV ") + latb + " ADD");
+  fw.SetParameter("lonlat2y", MString("(%lat-") + latb + ")*60");
+  fw.SetParameter("lonlatuv2dotx", "%u*(0.864/1.852)/Cos(%lat*M_PI/180)");
+  fw.SetParameter("lonlatuv2doty", "%v*(0.864/1.852)");
+  fw.SetParameter("dotxdotylonlat2u_rp", "EXCH POP EXCH POP COSD 1.852 0.864 DIV MUL MUL");
+  fw.SetParameter("dotxdotylonlat2v_rp", "POP POP 1.852 0.864 DIV MUL EXCH POP");
+  fw.SetParameter("dotxdotylonlat2u", "(1.852/0.864)*Cos(%lat*M_PI/180)*%dotx");
+  fw.SetParameter("dotxdotylonlat2v", "(1.852/0.864)*%doty");
+  fw.SetParameter("BeginDate", times.front().ToTString());
+  fw.SetParameter("EndDate", times.back().ToTString());
+  fw.SetParameter("Timestep", (times[1] - times[0]).Seconds());
+  fw.SetParameter("DataID", michlib::UniqueId());
+  fw.SetParameter("Creation command", args.at("_cmdline"));
+  fw.SetParameter("U variable", uvar);
+  fw.SetParameter("V variable", vvar);
+ }
+
+ for(size_t it = 0; it < times.size(); it++)
+ {
+  michlib::message(times[it].ToTString());
+
+  const bool is2D = ad.Is2D(uvar);
+
+  std::shared_ptr<Data2D> udata2, vdata2;
+  std::shared_ptr<Data3D> udata3, vdata3;
+
+  auto U = [&data2 = std::as_const(udata2), &data3 = std::as_const(udata3), is2D = is2D](size_t ix, size_t iy) { return is2D ? (*data2)(ix, iy) : (*data3)(ix, iy, 0); };
+  auto V = [&data2 = std::as_const(vdata2), &data3 = std::as_const(vdata3), is2D = is2D](size_t ix, size_t iy) { return is2D ? (*data2)(ix, iy) : (*data3)(ix, iy, 0); };
+
+  auto Ufill = [&data2 = std::as_const(udata2), &data3 = std::as_const(udata3), is2D = is2D](size_t ix, size_t iy)
+  { return is2D ? data2->IsFill(ix, iy) : data3->IsFill(ix, iy, 0); };
+  auto Vfill = [&data2 = std::as_const(vdata2), &data3 = std::as_const(vdata3), is2D = is2D](size_t ix, size_t iy)
+  { return is2D ? data2->IsFill(ix, iy) : data3->IsFill(ix, iy, 0); };
+
+  auto IsCoast = [&U = Ufill, &V = Vfill](size_t ix, size_t iy) { return U(ix, iy) || V(ix, iy); };
+
+  if(is2D)
+  {
+   auto retu = ad.Read2D(uvar, ind[it]);
+   auto retv = ad.Read2D(vvar, ind[it]);
+   if(!retu) return "Can't read U data";
+   if(!retv) return "Can't read V data";
+   udata2 = retu.Value();
+   vdata2 = retv.Value();
+  }
+  else
+  {
+   auto retu = ad.Read3D(uvar, ind[it]);
+   auto retv = ad.Read3D(vvar, ind[it]);
+   if(!retu) return "Can't read U data";
+   if(!retv) return "Can't read V data";
+   udata3 = retu.Value();
+   vdata3 = retv.Value();
+  }
+
+  for(size_t ilat = 0; ilat < Ny; ilat++)
+  {
+   for(size_t ilon = 0; ilon < Nx; ilon++)
+   {
+    if(IsCoast(ilon, ilat))
+    {
+     fw.Write(fillval);
+     fw.Write(fillval);
+    }
+    else
+    {
+     fw.Write(OutU(U(ilon, ilat)) * (0.864 / 1.852) / Cos(proj->Lat(ilon, ilat) * M_PI / 180.0));
+     fw.Write(OutV(V(ilon, ilat)) * (0.864 / 1.852));
+    }
+   }
+  }
+ }
+
+ return "";
+}
+
+template<class D>
+requires(!AdapterSupported<D> && ReadIsGrid<D>)
+MString ActionGenIntFile_DoAction(const CLArgs& args, D& ds)
 {
  auto [reg, regerr] = GetRegion<D>(args);
  if(regerr.Exist()) return regerr;

actions/actioninfo.h

@@ -24,10 +24,15 @@ MString ActionInfo_DoAction(const CLArgs& args, D& data)
   ad = std::move(ret.Value());
  }
 
- { // Modules
+ // Modules
+ {
   auto ret = CF::TEOS10Module(args, pars, ad);
   if(!ret) return "Error adding TEOS-10 derived variables";
  }
+ {
+  auto ret = CF::GeostrophicModule(args, pars, ad);
+  if(!ret) return "Error adding geostrophic velocity";
+ }
 
  const MString titCol   = F(C::BWHITE);
  const MString timeCol  = F(C::YELLOW);
@@ -56,6 +61,21 @@ MString ActionInfo_DoAction(const CLArgs& args, D& data)
   if(uniform) out += ", step is " + timeCol + (times.back() - times.front()).Seconds() / (times.size() - 1) + C() + " s";
   message(out);
  }
+ {
+  const auto times = ad.TimeIndexes() | std::views::transform([&times = ad.Times()](size_t i) { return times[i]; });
+
+  const bool onetime = times.size() == 1;
+  const bool uniform = CF::CheckUniform(times);
+
+  MString out = MString("Selected time moments: ") + timeCol + times.size() + C();
+  if(onetime)
+   out += ", " + timeCol + times[0].ToTString() + C();
+  else
+   out += ", from " + timeCol + times.front().ToTString() + C() + " to " + timeCol + times.back().ToTString() + C();
+  if(!onetime && !uniform) out += ", not uniform";
+  if(uniform) out += ", step is " + timeCol + (times.back() - times.front()).Seconds() / (times.size() - 1) + C() + " s";
+  message(out);
+ }
 
  // Variables
  std::map<std::shared_ptr<Projection>, MString> usedproj;
@@ -72,7 +92,7 @@ MString ActionInfo_DoAction(const CLArgs& args, D& data)
   // Metainfo
   message(ad.Is2D(name) ? v2dCol : v3dCol, name, C(), info->lname.Exist() ? (" - " + info->lname) : "",
           info->targetunit.Exist() ? (", measured in " + unitCol + info->targetunit + C()) : "");
-  if(info->comment.Exist()) message(comCol," ", info->comment, C());
+  if(info->comment.Exist()) message(comCol, " ", info->comment, C());
   if(info->stname.Exist()) message(" Standard name: " + info->stname);
 
   // Vertical grid description

actions/actiontsc.h

@@ -29,10 +29,15 @@ MString ActionTSC_DoAction(const CLArgs& args, D& ds)
   ad = std::move(ret.Value());
  }
 
- { // Modules
+ // Modules
+ {
   auto ret = CF::TEOS10Module(args, pars, ad);
   if(!ret) return "Error adding TEOS-10 derived variables";
  }
+ {
+  auto ret = CF::GeostrophicModule(args, pars, ad);
+  if(!ret) return "Error adding geostrophic velocity";
+ }
 
  MString name      = args.contains("out") ? args.at("out") : "out.nc";
  int     compress  = args.contains("compress") ? args.at("compress").ToInt() : 3;
@@ -65,9 +70,6 @@ MString ActionTSC_DoAction(const CLArgs& args, D& ds)
   }
  }
 
- // Delete unused variables
- ad.CleanVariables(vars);
-
  pars.SetParameter("variables", varstring);
 
  NCFileW ncfw;
@@ -75,7 +77,7 @@ MString ActionTSC_DoAction(const CLArgs& args, D& ds)
  {
   decltype(pars) emptypars;
   const auto&    upars = obfuscate ? emptypars : pars;
-  auto           ret   = ncfw.Create(name, ad, upars, compress);
+  auto           ret   = ncfw.Create(name, ad, vars, upars, compress);
   if(!ret) return "Can't create file " + name;
  }
 
@@ -85,14 +87,14 @@ MString ActionTSC_DoAction(const CLArgs& args, D& ds)
  {
   michlib::message(ad.Times()[ind[it]].ToTString());
   for(const auto& v: ad.Vars())
-   if(ad.Is2D(v.first))
+   if(vars.contains(v.first) && ad.Is2D(v.first))
    {
     auto data = ad.Read2D(v.first, ind[it]);
     if(!data) return "Can't read data";
     auto ret = ncfw.WriteVariable(v.first, *data.Value(), it);
     if(!ret) return "Can't write data";
    }
-   else
+   else if(vars.contains(v.first) && ad.Is3D(v.first))
    {
     auto data = ad.Read3D(v.first, ind[it]);
     if(!data) return "Can't read data";

doc.txt

@@ -34,7 +34,8 @@ int - интерполирует данные для точек траектор
 genintfile - делает файл со скоростями для расчёта лагранжевых траекторий.
  time, timeb, timee - смотри описание tsc.
  lonb, lone, latb, late - смотри описание tsc.
- geostrophic. Смотри описание этого параметра для action=uv.
+ geostrophic. Смотри описание этого параметра для action=uv. Игнорируется для источников HYCOM и COPERNICUS.
+ u, v - имена переменных для зональной и меридиональной скоростей. По умолчанию - auto (выбираются имена "u", "v" или, при отсутствии таких переменных "ugs", "vgs"). Работает только для источников HYCOM и COPERNICUS.
  out - выходной файл.
 grad - Работает только с источником TSCDATA. Берётся выходной файл от action=tsc и для всех переменных в нём рассчитывает градиенты по методу Белкина.
  out - выходной файл.
@@ -52,8 +53,8 @@ BINFILE - файл со скоростями для расчёта лагран
 MODISBINLOCAL - локальное зеркало данных по температуре поверхности моря и концентрации хлорофилла спутников Aqua- и Terra-MODIS. Нужно указывать регион.
  dataset - одно или несколько ключевых слов, разделённых запятыми, из списка A_CHL, T_CHL, CHL (A_CHL, T_CHL), A_SST, T_SST, A_SST4, T_SST4, A_NSST, T_NSST, A_DSST (A_SST, A_SST4), T_DSST (T_SST, T_SST4), A_ALL (A_DSST, A_NSST), T_ALL (T_DSST, T_NSST), SST (A_SST, T_SST), SST4 (A_SST4, T_SST4), DSST (A_DSST, T_DSST), NSST (A_NSST, T_NSST), ALL (A_ALL, T_ALL). Ключевые слова, соответствующие температуре и хлорофиллу, не должны смешиваться и должны соответствовать запрашиваемой переменной (chl или temp). По умолчанию CHL или ALL, в зависимости от параметра var.
 HYCOM - данные реанализа по модели HYCOM. Нужно указывать регион.
- dataset - Forecast или Hindcast. По умолчанию Forecast.
- layer или depth - выбор слоя или глубины. Подробнее смотри описание этих параметров для источника NEMO.
+ dataset - набор данных, задаваемый в файле конфигурации.
+ Регион, слои, время и переменные задаются так же, как и для COPERNICUS.
 NEMO - данные реанализа по модели NEMO. Нужно указывать регион.
  dataset - DT, NRT, NRT6 (глобальные), BALTICDT, BALTICNRT, BALTICNRT1, BLKSEADT, BLKSEANRT, MEDSEADT, MEDSEANRT, BISCDT, BISCNRT, ENWSDT, ENWSNRT (региональные). По умолчанию DT.
  layer или depth - выбор номера слоя или глубины в метрах. Если указаны оба параметра, приоритет имеет depth. Если указан depth, выбирается слой, ближайший к заданной глубине. Нумерация слоёв начинается с нуля и идёт сверху вниз. Если не указаны оба параметра, используется нулевой слой (поверхность). Эти параметры игнорируются для некоторых переменных, например, ssh или mld.

include/Adapter.h

@@ -113,11 +113,6 @@ class Adapter
   return true;
  }
 
- void CleanVariables(const std::set<MString>& keep)
- {
-  std::erase_if(vars, [&keep](const auto& k) { return !keep.contains(k.first); });
- }
-
  const auto& Times() const { return times; }
  const auto& TimeIndexes() const { return tindexes; }
  const auto& Vars() const { return vars; }

include/CF.h

@@ -28,6 +28,19 @@ class CF
   virtual ~VInfoTEOS10() override = default;
  };
 
+ struct VInfoGeostrophic: public Adapter::VInfoBase
+ {
+  enum Operation
+  {
+   U,
+   V
+  };
+  MString   sshvar;
+  Operation op;
+
+  virtual ~VInfoGeostrophic() override = default;
+ };
+
  // Get time variable name
  static MString GetTVarName(const NCZarr& nc);
 
@@ -94,6 +107,8 @@ class CF
  }
 
  static RetVal<std::shared_ptr<Data3D>> TEOS103DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it);
+ static RetVal<std::shared_ptr<Data2D>> Geostrophic2DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr<Projection> proj, size_t it);
 
  static Error TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad);
+ static Error GeostrophicModule(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad);
 };

include/ncfilew.h

@@ -389,10 +389,10 @@ class NCFileW: public NCFileWBase
 
  michlib::Error AddTimeData(const std::vector<MDateTime>& times);
 
- michlib::Error Create(const MString& name, const Adapter& ad, const michlib_internal::ParameterListEx& pars, int compression);
+ michlib::Error Create(const MString& name, const Adapter& ad, const std::set<MString>& vars, const michlib_internal::ParameterListEx& pars, int compression);
 
  template<class D>
- requires (std::is_same_v<D, Data2D> || std::is_same_v<D, Data3D>)
+ requires(std::is_same_v<D, Data2D> || std::is_same_v<D, Data3D>)
  michlib::Error WriteVariable(const MString& name, const D& data, size_t it);
 
  template<class D> MString Create(const D& data, const MString& name, int compression)

sources/HYCOM.cpp

@@ -0,0 +1,117 @@
+#define MICHLIB_NOSOURCE
+#include "HYCOM.h"
+#include "CF.h"
+#include <numeric>
+#include <valarray>
+
+using michlib::GPL;
+
+RetVal<Adapter> HYCOMData::GetAdapter(const CLArgs& args, michlib_internal::ParameterListEx& pars) const
+{
+ static const MString pref = "HYCOMData::GetAdapter";
+
+ //const bool debug = args.contains("debug");
+
+ Adapter out(Adapter::ZarrMethod::MNC);
+
+ GPL.UsePrefix("HYCOM");
+
+ // Get link for dataset
+ const MString dataset = args.contains("dataset") ? args.at("dataset") : "Forecast";
+
+ pars.AddParameter("dataset", dataset);
+
+ std::vector<MString> urls;
+ size_t               cachettl = GPL.ParameterUValue(dataset + "_CacheTTL", 3600 * 24 * 7);
+ {
+  const auto url = GPL.ParameterSValue(dataset + "_URL", "");
+  if(url.Exist())
+   urls.push_back(url);
+  else
+  {
+   size_t i = 1;
+   while(true)
+   {
+    const auto url = GPL.ParameterSValue(dataset + "_URL" + i, "");
+    if(!url.Exist()) break;
+    urls.push_back(url);
+    i++;
+   }
+  }
+ }
+ if(urls.size() == 0) return {pref, "No urls find for datest " + dataset};
+
+ out.SetTitle(GPL.ParameterSValue(dataset + "_Title", "Unknown"));
+ pars.AddParameter("title", out.Title());
+
+ // Open data
+ std::unique_ptr<NCZarr> pnc;
+ std::vector<MDateTime>  times;
+ {
+  std::vector<Adapter::TimeRow> table;
+  for(const auto& url: urls)
+  {
+   Adapter::TimeRow       row;
+   std::vector<MDateTime> ltimes;
+
+   row.beg         = table.size() == 0 ? 0 : (table.back().end + 1);
+   row.file        = url;
+   auto [val, res] = cache->Get(url);
+   if(!res)
+   {
+    pnc.reset(new NCZarr);
+    auto ret = pnc->OpenMultiNC(url.Split(";"));
+    if(!ret) return ret.Add(pref, "Can't open link " + url);
+    auto tvar = CF::GetTVarName(*pnc);
+    if(!tvar.Exist()) return {pref, "Can't find time variable in the dataset " + dataset};
+    if(pnc->NDim(tvar) != 1) return {pref, "Unsupported number of time dimensions: " + MString(pnc->NDim(tvar))};
+
+    auto timesret = CF::ReadTimes(*pnc, tvar);
+    if(!timesret) return ret.Add(pref, "Can't read time values from the dataset " + dataset);
+    ltimes = timesret.Value();
+    if(ltimes.size() == 0) return {pref, "No time moments in " + url};
+
+    michlib::Serializer s;
+    s.Serialize(ltimes.size());
+    for(size_t i = 0; i < ltimes.size(); i++) s.Serialize(ltimes[i]);
+    val = MString(s.Buf(), s.Len());
+    cache->Put(url, val, cachettl);
+   }
+   else
+   {
+    michlib::Deserializer ds(val.Buf(), val.Len());
+    size_t                vsize;
+    bool                  ok;
+    ok = ds.Deserialize(vsize);
+    if(!ok) return {pref, "Cache line for url " + url + " is bad"};
+    ltimes.resize(vsize);
+    for(size_t i = 0; i < vsize; i++) ok = ok && ds.Deserialize(ltimes[i]);
+    if(!ok) return {pref, "Cache line for url " + url + " is bad"};
+   }
+   row.end = row.beg + ltimes.size() - 1;
+   times.insert(times.end(), ltimes.begin(), ltimes.end());
+   table.push_back(row);
+  }
+  out.SetTimeTable(std::move(table));
+ }
+
+ if(!pnc)
+ {
+  pnc.reset(new NCZarr);
+  auto ret = pnc->OpenMultiNC(urls[0].Split(";"));
+  if(!ret) return ret.Add(pref, "Can't open link " + urls[0]);
+ }
+
+ // Set times
+ {
+  auto ret = out.SetTimes(std::move(times), args, pars);
+  if(!ret) return ret.Add(pref, "Can't filter time values from the dataset " + dataset);
+ }
+
+ {
+  auto ret = CF::FillAdapterFromCF(args, pars, out, std::move(pnc));
+  if(!ret) return ret.Add(pref, "Can't get adapter");
+ }
+
+ return out;
+}

sources/HYCOM.h

@@ -14,10 +14,23 @@ class HYCOMData: public LayeredData
 
  Type type = TYPE_UNKNOWN;
 
+ std::unique_ptr<GenericCache> cache;
+
  public:
  static constexpr const char* name = "HYCOM";
 
- HYCOMData() = default;
+ HYCOMData()
+ {
+  auto oldprefix = michlib::GPL.UsePrefix("HYCOM");
+  // Cache
+  cache.reset(CreateCache(michlib::GPL.ParameterSValue("Cache", "")));
+  if(!cache)
+  {
+   michlib::errmessage("Can't init cache");
+   cache.reset(new FakeCache);
+  }
+  michlib::GPL.UsePrefix(oldprefix);
+ }
 
  MString DataTitle() const
  {
@@ -51,4 +64,7 @@ class HYCOMData: public LayeredData
   SetTitle(DataTitle());
   return LayeredData::Open(dataset);
  }
+
+ // Adapter
+ RetVal<Adapter> GetAdapter(const CLArgs& args, michlib_internal::ParameterListEx& pars) const;
 };

src/Adapter.cpp

@@ -256,7 +256,7 @@ Error Adapter::ReadCheck(const MString& var, size_t it)
   }
   case(ZarrMethod::MNC):
   {
-   auto ret = nc->OpenMultiNC(curfile.Split(":"));
+   auto ret = nc->OpenMultiNC(curfile.Split(";"));
    if(!ret) return ret.Add(pref, "OpenMultiNC failed for file " + curfile);
    break;
   }

src/CF.cpp

@@ -117,7 +117,7 @@ RetVal<std::vector<MDateTime>> CF::ReadTimes(const NCZarr& nc, const MString& tn
  }
 
  out.resize(time.size());
- for(size_t i = 0; i < time.size(); i++) out[i] = refdate + static_cast<time_t>(time[i] * step);
+ for(size_t i = 0; i < time.size(); i++) out[i] = refdate + static_cast<time_t>(michlib::Round(time[i])) * step;
  return out;
 }
 
@@ -256,7 +256,10 @@ Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListE
   }
 
   // Check uniform and create projection
-  if(CF::CheckUniform(rlons, 0.0005) && CF::CheckUniform(rlats, 0.0005))
+  bool lonuni = CF::CheckUniform(rlons, 0.005);
+  bool latuni = CF::CheckUniform(rlats, 0.005);
+
+  if(lonuni && latuni)
   {
    proj.reset(new Projection(Projection::Create<Projection::Type::EQC>(rlons, rlats)));
    pars.AddParameter("lonstep", (rlons.back() - rlons.front()) / (rlons.size() - 1));
@@ -264,7 +267,8 @@ Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListE
   }
   else
   {
-   michlib::errmessage("Warning! Nonuniform grid");
+   if(!lonuni) michlib::errmessage("Warning! Nonuniform longitude grid");
+   if(!latuni) michlib::errmessage("Warning! Nonuniform latitude grid");
    proj.reset(new Projection(Projection::Create<Projection::Type::IRREGULARXY>(std::move(rlons), std::move(rlats))));
   }
  }
@@ -409,12 +413,13 @@ Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListE
  // Read variables information
  for(const auto& v: pnc->VarNames())
  {
-  auto vinfo  = new Adapter::VInfoNC();
   auto dnames = pnc->DimNames(v);
   if(dnames.size() != 3 && dnames.size() != 4) continue;
   if(dnames.size() == 3 && (dnames[0] != axisinfo.tdimname || dnames[1] != axisinfo.ydimname || dnames[2] != axisinfo.xdimname)) continue;
   if(dnames.size() == 4 && (dnames[0] != axisinfo.tdimname || dnames[1] != axisinfo.zdimname || dnames[2] != axisinfo.ydimname || dnames[3] != axisinfo.xdimname)) continue;
 
+  auto vinfo = new Adapter::VInfoNC();
+
   vinfo->name = v;
   vinfo->unit = "1";
   if(dnames.size() == 4)
@@ -437,6 +442,7 @@ Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListE
       },
       pnc->VarFill(v));
 
+  if(vinfo->unit == "psu") vinfo->unit = "1e-3"; // Deal with salinity psu
   vinfo->targetunit = vinfo->unit;
   if(vinfo->stname.SubStr(vinfo->stname.Len() - 7, 8) == "velocity") vinfo->targetunit = "cm/s";
 
@@ -450,7 +456,7 @@ Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListE
 
 RetVal<std::shared_ptr<Data3D>> CF::TEOS103DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it)
 {
- static const MString pref = "CF::TempSalDer3DReader";
+ static const MString pref = "CF::Teos103DReader";
 
  const auto* v = dynamic_cast<const struct VInfoTEOS10*>(vinfo);
  if(v == nullptr) return {pref, "Incorrect argument type"};
@@ -544,9 +550,119 @@ RetVal<std::shared_ptr<Data3D>> CF::TEOS103DReader(Adapter& ad, const Adapter::V
  return out;
 }
 
+RetVal<std::shared_ptr<Data2D>> CF::Geostrophic2DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr<Projection> proj, size_t it)
+{
+ static const MString pref = "CF::Geostrophic2DReader";
+
+ const auto* v = dynamic_cast<const struct VInfoGeostrophic*>(vinfo);
+ if(v == nullptr) return {pref, "Incorrect argument type"};
+
+ std::shared_ptr<Data2D> out(new Data2D(proj, {.unit = v->targetunit, .stname = v->stname, .lname = v->lname, .comment = v->comment}));
+
+ std::shared_ptr<Data2D> ssh;
+ {
+  auto retssh = ad.Read2D(v->sshvar, ad.CurrentTime());
+  if(!retssh) return {pref, "Can't read ssh variable " + v->sshvar};
+  ssh = retssh.Value();
+ }
+
+ if(ssh->Ny() != out->Ny()) return {pref, MString("For ssh variable Ny=") + ssh->Ny() + ", for output variable Ny=" + out->Ny()};
+ if(ssh->Nx() != out->Nx()) return {pref, MString("For ssh variable Nx=") + ssh->Nx() + ", for output variable Nx=" + out->Nx()};
+
+ const UtUnit infr(ssh->Unit()), into("cm");
+ const UtUnit outfr("cm/s"), outto(out->Unit());
+ const auto   cnvin  = Converter(infr, into);
+ const auto   cnvout = Converter(outfr, outto);
+
+ if(!cnvin) return {pref, "Can't convert units from \"" + ssh->Unit() + "\" to \"cm\""};
+ if(!cnvout) return {pref, "Can't convert units from \"cm/s\" to \"" + ssh->Unit() + "\""};
+
+ const bool needsconvertin  = infr != into;
+ const bool needsconvertout = outfr != outto;
+
+ auto In  = [needsconvertin, &cnvin = cnvin](real r) -> real { return needsconvertin ? cnvin(r) : r; };
+ auto Out = [needsconvertout, &cnvout = cnvout](real r) -> real { return needsconvertout ? cnvout(r) : r; };
+
+ switch(v->op)
+ {
+ case(VInfoGeostrophic::U):
+ {
+  for(size_t iy = 0; iy < ssh->Ny(); iy++)
+  {
+   const real phi   = michlib::M_PI * ssh->Lat(iy) / 180.0;
+   const real sphi  = michlib::Sin(phi);
+   const real s2phi = michlib::Sin(2.0 * phi);
+   const real g     = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi);                   // g in cm/s
+   const real f     = 2.0 * 7.29211 * sphi;                                                                   // Coriolis parameter multiplied by 100'000
+   const real db    = (iy == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh->Lat(iy) - ssh->Lat(iy - 1)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter
+   const real du =
+       (iy == ssh->Ny() - 1) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh->Lat(iy + 1) - ssh->Lat(iy)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter
+
+   for(size_t ix = 0; ix < ssh->Nx(); ix++)
+   {
+    const bool bok = iy != 0 && !ssh->IsFill(ix, iy - 1) && !ssh->IsFill(ix, iy);
+    const bool uok = iy != ssh->Ny() - 1 && !ssh->IsFill(ix, iy + 1) && !ssh->IsFill(ix, iy);
+    if(!(bok || uok))
+     (*out)(ix, iy) = out->Fillval();
+    else if(!bok)
+     (*out)(ix, iy) = Out(-g * (In((*ssh)(ix, iy + 1)) - In((*ssh)(ix, iy))) / (f * du));
+    else if(!uok)
+     (*out)(ix, iy) = Out(-g * (In((*ssh)(ix, iy)) - In((*ssh)(ix, iy - 1))) / (f * db));
+    else
+    {
+     const real bder = -g * (In((*ssh)(ix, iy)) - In((*ssh)(ix, iy - 1))) / (f * db);
+     const real uder = -g * (In((*ssh)(ix, iy + 1)) - In((*ssh)(ix, iy))) / (f * du);
+     (*out)(ix, iy)  = Out(0.5 * (bder + uder));
+    }
+   }
+  }
+  break;
+ }
+ case(VInfoGeostrophic::V):
+ {
+  for(size_t iy = 0; iy < ssh->Ny(); iy++)
+  {
+   const real phi   = michlib::M_PI * ssh->Lat(iy) / 180.0;
+   const real sphi  = michlib::Sin(phi);
+   const real cphi  = michlib::Cos(phi);
+   const real s2phi = michlib::Sin(2.0 * phi);
+   const real g     = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi); // g in cm/s
+   const real f     = 2.0 * 7.29211 * sphi;                                                 // Coriolis parameter multiplied by 100'000
+
+   for(size_t ix = 0; ix < ssh->Nx(); ix++)
+   {
+    const real dl =
+        (ix == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh->Lon(ix) - ssh->Lon(ix - 1)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter
+    const real dr = (ix == ssh->Nx() - 1)
+                        ? 0.0
+                        : (michlib::M_PI * 6371.0 * (ssh->Lon(ix + 1) - ssh->Lon(ix)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter
+
+    const bool lok = ix != 0 && !ssh->IsFill(ix - 1, iy) && !ssh->IsFill(ix, iy);
+    const bool rok = ix != ssh->Nx() - 1 && !ssh->IsFill(ix + 1, iy) && !ssh->IsFill(ix, iy);
+    if(!(lok || rok))
+     (*out)(ix, iy) = out->Fillval();
+    else if(!lok)
+     (*out)(ix, iy) = Out(g * (In((*ssh)(ix + 1, iy)) - In((*ssh)(ix, iy))) / (f * dr));
+    else if(!rok)
+     (*out)(ix, iy) = Out(g * (In((*ssh)(ix, iy)) - In((*ssh)(ix - 1, iy))) / (f * dl));
+    else
+    {
+     const real lder = g * (In((*ssh)(ix, iy)) - In((*ssh)(ix - 1, iy))) / (f * dl);
+     const real rder = g * (In((*ssh)(ix + 1, iy)) - In((*ssh)(ix, iy))) / (f * dr);
+     (*out)(ix, iy)  = Out(0.5 * (lder + rder));
+    }
+   }
+  }
+  break;
+ }
+ }
+
+ return out;
+}
+
 Error CF::TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad)
 {
- static const MString pref = "TempSalDerModule";
+ static const MString pref = "TEOS10Module";
 
  MString sal, temp, ptemp;
 
@@ -702,3 +818,62 @@ Error CF::TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pa
 
  return {};
 }
+
+Error CF::GeostrophicModule(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad)
+{
+ static const MString pref = "GeostrophicModule";
+
+ MString ssh;
+
+ for(const auto& v: ad.Vars())
+ {
+  if(!ad.Is2D(v.first)) continue;
+  const auto& info = v.second.info;
+  if(!ssh.Exist() && info->stname == "sea_surface_height_above_geoid") ssh = v.first;
+  if(!ssh.Exist() && info->stname == "sea_surface_elevation") ssh = v.first;
+ }
+
+ if(!ssh.Exist()) return {};
+
+ const auto& svar  = ad.Vars().at(ssh);
+ const auto& sproj = svar.proj;
+ const auto& rinfo = svar.info->rinfo;
+
+ if(!sproj->IsRectangular()) return {}; // We deal only with rectangular grids for now
+
+ { // U-component
+  auto vinfo = new VInfoGeostrophic();
+
+  vinfo->sshvar = ssh;
+  vinfo->op     = VInfoGeostrophic::U;
+
+  vinfo->stname     = "surface_geostrophic_eastward_sea_water_velocity";
+  vinfo->lname      = "Geostrophic eastward velocity";
+  vinfo->comment    = "Module Geostrophic: Geostrophic eastward velocity calculated from sea surface height (" + ssh + ")";
+  vinfo->targetunit = "cm/s";
+  vinfo->rinfo      = rinfo;
+  vinfo->fillval    = NAN;
+
+  const bool addsuc = ad.Add2DVariable("ugs", Adapter::VInfo(vinfo), sproj, Geostrophic2DReader);
+  if(!addsuc) return {pref, "Variable with name ugs already exist"};
+ }
+
+ { // V-component
+  auto vinfo = new VInfoGeostrophic();
+
+  vinfo->sshvar = ssh;
+  vinfo->op     = VInfoGeostrophic::V;
+
+  vinfo->stname     = "surface_geostrophic_northward_sea_water_velocity";
+  vinfo->lname      = "Geostrophic northward velocity";
+  vinfo->comment    = "Module Geostrophic: Geostrophic northward velocity calculated from sea surface height (" + ssh + ")";
+  vinfo->targetunit = "cm/s";
+  vinfo->rinfo      = rinfo;
+  vinfo->fillval    = NAN;
+
+  const bool addsuc = ad.Add2DVariable("vgs", Adapter::VInfo(vinfo), sproj, Geostrophic2DReader);
+  if(!addsuc) return {pref, "Variable with name vgs already exist"};
+ }
+
+ return {};
+}

src/copcat.cpp

@@ -1,6 +1,7 @@
 #define MICHLIB_NOSOURCE
 #include "copcat.h"
 #include "GPL.h"
+#include "conscolors.h"
 #include "mirrorfuncs.h"
 
 const MString CopernicusCatalog::caturl = "https://stac.marine.copernicus.eu/metadata/catalog.stac.json";
@@ -205,7 +206,7 @@ RetVal<Json::Value> CopernicusCatalog::GetJSON(const MString& url) const
   content = std::move(val);
  else
  {
-  michlib::message(url + " not found in cache, downloading");
+  michlib::message(F(C::BBLACK), url, " not found in cache, downloading", C());
   auto [out, res] = GetUrl(chandle, url);
   if(res != CURLE_OK) return Error(pref, MString("can't download JSON: ") + chandle.Err());
   cache->Put(url, out, 3600);

src/ncfilew.cpp

@@ -184,7 +184,7 @@ michlib::Error NCFileW::AddTimeData(const std::vector<MDateTime>& times)
  return {};
 }
 
-michlib::Error NCFileW::Create(const MString& name, const Adapter& ad, const michlib_internal::ParameterListEx& pars, int compression)
+michlib::Error NCFileW::Create(const MString& name, const Adapter& ad, const std::set<MString>& vars, const michlib_internal::ParameterListEx& pars, int compression)
 {
  static const MString pref = "NCFileW::Create";
  if(Ok()) return {pref, "File already opened"};
@@ -211,6 +211,7 @@ michlib::Error NCFileW::Create(const MString& name, const Adapter& ad, const mic
 
  for(const auto& v: ad.Vars())
  {
+  if(!vars.contains(v.first)) continue;
   bool is2D = v.second.Read2D ? true : false;
   bool is3D = v.second.Read3D ? true : false;
 

src/zarr.cpp

@@ -1,5 +1,6 @@
 #define MICHLIB_NOSOURCE
 #include "zarr.h"
+#include "conscolors.h"
 #include "copcat.h"
 #include <blosc.h>
 
@@ -171,7 +172,7 @@ Error ZarrFunctions::GetChunk(const MString& var, const std::vector<size_t>& chu
 
  if(!suc)
  {
-  michlib::message(str + " not found in cache, downloading");
+  michlib::message(F(C::BBLACK), str, " not found in cache, downloading", C());
   CURLRAII myhandle; // TODO: remove this workaround of unknown bug
   if(proxyurl.Exist()) curl_easy_setopt(myhandle, CURLOPT_PROXY, proxyurl.Buf());
   //auto [out, res] = GetUrl(chandle, str);