Generate usage information from doc.txt

CMakeLists.txt

@@ -34,7 +34,7 @@ set(BoldWhite   "${Esc}[1;37m")
 set(STATIC_ANALYZER OFF CACHE BOOL "Using GCC static analyzer if available")
 
 project(odm C CXX)
-set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD 23)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_INTERPROCEDURAL_OPTIMIZATION ON)
 

GSW-C

@@ -1 +1 @@
-Subproject commit 83c1eb7503a6eac8231f3cde4e10f17aed52f30f
+Subproject commit e43a95cd5c8aad39895d8f57fdbdcc9aced8657e

actions/actioninfo.h

@@ -1,12 +1,154 @@
 #pragma once
+#include "CF.h"
 #include "actiondep.h"
-#include "merrors.h"
 
 using michlib::message;
+using michlib::uint2;
 
-ADD_ACTION(Info, info, InfoSupported<Source>);
+ADD_ACTION(Info, info, (InfoSupported<Source> || AdapterSupported<Source>));
 
-template<class D> MString ActionInfo::DoAction(const CLArgs& args, D& data)
+template<class D> MString ActionInfo_DoAction(const CLArgs& args, D& data);
+
+template<class D> MString ActionInfo::DoAction(const CLArgs& args, D& data) { return ActionInfo_DoAction(args, data); }
+
+template<class D>
+requires AdapterSupported<D>
+MString ActionInfo_DoAction(const CLArgs& args, D& data)
+{
+ Adapter                           ad;
+ michlib_internal::ParameterListEx pars;
+ pars.UsePrefix("");
+ {
+  auto ret = data.GetAdapter(args, pars);
+  if(!ret) return "Error opening source";
+  ad = std::move(ret.Value());
+ }
+
+ const MString titCol   = F(C::BWHITE);
+ const MString timeCol  = F(C::YELLOW);
+ const MString v2dCol   = F(C::BGREEN);
+ const MString v3dCol   = F(C::CYAN);
+ const MString unitCol  = F(C::YELLOW);
+ const MString usedCol  = F(C::BMAGENTA);
+ const MString depthCol = F(C::YELLOW);
+ const MString horCol   = F(C::YELLOW);
+
+ // Header
+ message(titCol, ad.Title(), C());
+ {
+  const auto& times = ad.Times();
+
+  const bool onetime = times.size() == 1;
+  const bool uniform = CF::CheckUniform(times);
+
+  MString out = MString("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;
+ std::map<std::shared_ptr<Vertical>, MString>   usedvert;
+ for(const auto& v: ad.Vars())
+ {
+  const auto& name = v.first;
+  const auto& info = v.second.info;
+  const auto& proj = v.second.proj;
+  const auto& vert = v.second.vert;
+
+  message("");
+
+  // 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.stname.Exist()) message(" Standard name: " + info.stname);
+
+  // Vertical grid description
+  if(!ad.Is2D(name))
+  {
+   if(usedvert.contains(vert))
+    message(" Vertical grid is same as for ", usedCol, usedvert[vert], C());
+   else
+   {
+    MString out = " Vertical grid: ";
+    switch(vert->Type())
+    {
+    case(Vertical::Type::LAYERS):
+    {
+     out += depthCol + vert->Nz() + C() + " layers";
+     break;
+    }
+    case(Vertical::Type::HORIZONTS):
+    {
+     out += depthCol + vert->Nz() + C() + " horizonts,";
+     for(size_t i = 0; i < vert->Nz(); i++) { out += " (" + MString(i) + " " + depthCol + vert->Depth(i) + C() + " m)"; }
+     break;
+    }
+    case(Vertical::Type::GRID):
+    {
+     out += depthCol + vert->Nz() + C() + " horizonts, from " + depthCol + vert->Depth(0) + C() + "m to " + depthCol + vert->Depth(vert->Nz() - 1) + C() + "m, step " + depthCol +
+            (vert->Depth(vert->Nz() - 1) - vert->Depth(0)) / (vert->Nz() - 1) + C() + "m";
+     break;
+    }
+    case(Vertical::Type::FIXED):
+    {
+     out += depthCol + vert->Depth(0) + C() + " m";
+     break;
+    }
+    case(Vertical::Type::UNDEFINED):
+    {
+     out += "undefined";
+     break;
+    }
+    }
+    message(out);
+    usedvert[vert] = name;
+   }
+  }
+
+  // Horizontal grid description
+  if(usedproj.contains(proj))
+   message(" Horizontal grid is same as for ", usedCol, usedproj[proj], C());
+  else
+  {
+   MString out = " Horizontal grid: ";
+
+   switch(proj->Type())
+   {
+   case(Projection::Type::IRREGULAR):
+   {
+    out += "irregular, " + horCol + proj->N() + C() + " points";
+    break;
+   }
+   case(Projection::Type::IRREGULARXY):
+   {
+    out += "irregular, " + horCol + proj->Nx() + C() + " x " + horCol + proj->Ny() + C() + " points, region (" + horCol + proj->Lon(0) + C() + " : " + horCol +
+           proj->Lon(proj->Nx() - 1) + C() + ") x (" + horCol + proj->Lat(0) + C() + " : " + horCol + proj->Lat(proj->Ny() - 1) + C() + ")";
+    break;
+   }
+   case(Projection::Type::EQC):
+   {
+    out += "rectangular, " + horCol + proj->Nx() + C() + " x " + horCol + proj->Ny() + C() + " points (" + horCol + (proj->Lon(proj->Nx() - 1) - proj->Lon(0)) / (proj->Nx() - 1) +
+           C() + " x " + horCol + (proj->Lat(proj->Ny() - 1) - proj->Lat(0)) / (proj->Ny() - 1) + C() + "), region (" + horCol + proj->Lon(0) + C() + " : " + horCol +
+           proj->Lon(proj->Nx() - 1) + C() + ") x (" + horCol + proj->Lat(0) + C() + " : " + horCol + proj->Lat(proj->Ny() - 1) + C() + ")";
+    break;
+   }
+   }
+   message(out);
+   usedproj[proj] = name;
+  }
+ }
+ return "";
+}
+
+template<class D>
+requires(!AdapterSupported<D> && InfoSupported<D>)
+MString ActionInfo_DoAction(const CLArgs& args, D& data)
 {
  auto resop = data.Open(args);
  if(resop.Exist()) return "Can't open source: " + resop;

actions/actiontsc.h

@@ -7,9 +7,101 @@
 
 using michlib::BFileW;
 
-ADD_ACTION(TSC, tsc, ReadPSupported<Source> || ReadSupported<Source>);
+ADD_ACTION(TSC, tsc, AdapterSupported<Source> || ReadPSupported<Source> || ReadSupported<Source>);
 
-template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
+template<class D> MString ActionTSC_DoAction(const CLArgs& args, D& data);
+
+template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& data) { return ActionTSC_DoAction(args, data); }
+
+template<class D>
+requires AdapterSupported<D>
+MString ActionTSC_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());
+ }
+
+ MString name      = args.contains("out") ? args.at("out") : "out.nc";
+ int     compress  = args.contains("compress") ? args.at("compress").ToInt() : 3;
+ bool    obfuscate = args.contains("obfuscate");
+
+ MString varstring;
+ if(args.contains("var"))
+  varstring = args.at("var");
+ else
+ {
+  if(args.contains("vars"))
+   varstring = args.at("vars");
+  else
+  {
+   if constexpr(HasDefVars<D>)
+    varstring = ds.DefaultVars();
+   else
+    return "Variables not specified";
+  }
+ }
+
+ std::set<MString> vars;
+ {
+  auto vlist = michlib::Split_on_words(varstring, " ,", false);
+  for(const auto& vname: vlist)
+  {
+   if(vars.contains(vname)) return "Duplicate variable " + vname + " in list " + varstring;
+   if(!ad.Vars().contains(vname)) return "Variable " + vname + " does'nt supported by this dataset";
+   vars.insert(vname);
+  }
+ }
+
+ // Delete unused variables
+ ad.CleanVariables(vars);
+
+ pars.SetParameter("variables", varstring);
+
+ NCFileW ncfw;
+
+ {
+  decltype(pars) emptypars;
+  const auto&    upars = obfuscate ? emptypars : pars;
+  auto           ret   = ncfw.Create(name, ad, upars, compress);
+  if(!ret) return "Can't create file " + name;
+ }
+
+ const auto& ind = ad.TimeIndexes();
+
+ for(size_t it = 0; it < ind.size(); it++)
+ {
+  michlib::message(ad.Times()[ind[it]].ToTString());
+  for(const auto& v: ad.Vars())
+   if(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
+   {
+    auto data = ad.Read3D(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";
+   }
+ }
+
+ return "";
+}
+
+template<class D>
+requires(!AdapterSupported<D> && (ReadPSupported<D> || ReadSupported<D>))
+MString ActionTSC_DoAction(const CLArgs& args, D& ds)
 {
  auto [reg, regerr] = GetRegion<D>(args);
  if(regerr.Exist()) return regerr;

doc.txt

@@ -0,0 +1,85 @@
+color, nocolor - включить/выключить цветной вывод в консоль.
+
+Действия
+info - показывает информацию об источнике данных. Собственных параметров нет. Обычно нужен параметр dataset для указания источника данных.
+tsc - получает физические величины на некоторой плоскости из источника данных. Параметры:
+ var или vars - список разделённых запятыми физических величин. Возможные значения: U, U2, u, v, ugeo, vgeo, temp, ptemp, pdens, sal, chl, mld, ssh, w, NO3, PO4, Si, O2, prprod, Cchl. Доступность тех или иных величин зависит от источника данных. Если источник имеет список переменных по умолчанию, то этот параметр может быть опущен.
+ time или пара timeb, timee. Параметр time может быть как одним моментом времени, так и регулярным выражением. Также time, timeb и timee могут быть одним из слов BEGIN, BEG, FIRST (соответствуют минимальному времени, для которого доступны данные) или END, LAST (соответствуют максимальному времени, для которого доступны данные). Если time меньше, чем BEGIN, будет использоваться минимальное время, если time больше END - максимальное. Если time регулярное выражение, то из набора данных выбираются моменты времени соответствующие этому выражению, и по ним проводится усреднение. Пара timeb и timee работает по аналогичному принципу, только выбираются моменты времени между timeb и timee включительно.
+ lonb, lone, latb, late - регион, для которого берутся данные. Для некоторых источников данных обязателен, для других игнорируется, смотри описания источников.
+ out - выходной файл.
+ format - формат выходного файла. Может быть bin или nc (netcdf). Если этот параметр отсутствует, формат определяется по расширению выходного файла (.bin или .nc). Игнорируется для источника COPERNICUS.
+ compress - число от 0 до 9, степень сжатия netcdf файла. По умолчанию - 3.
+ gradient - если параметр присутствует, считать градиенты полей, а не сами поля. Игнорируется для источника COPERNICUS.
+ average - усреднить все поля по времени. Пока не реализовано для источника COPERNICUS.
+ obfuscate - если параметр присутствует, не записывать в выходные файлы метаинформацию.
+ debug - если параметр присутствует, записать в выходной файл дополнительную информацию. Игнорируется при наличии obfuscate.
+uv - получает скорость, производные от неё величины, вычисляет стационарные точки.
+ time, timeb, timee - смотри описание tsc.
+ lonb, lone, latb, late - смотри описание tsc.
+ geostrophic. Может быть быть пустым или любым из значений ssh, surf, surface, nongeo. Определяет, какую скорость из имеющихся в наборе данных использовать. Для значения ssh или пустого используются геострофические скорости, рассчитанные по полю высоты уровня моря. Для значений surf или surface используются гестрофические поверхностные скорости из данных. Для значения nongeo используются обычные скорости.
+ out - выходной файл со скоростями, дивергенцией, ротором, параметром Окуба-Вейсса, квадратом скорости и её дисперсии (вихревая кинетическая энергия). Если этот параметр отсутствует, соответствующий файл не будет записан.
+ outformat - формат выходного файла (out). Может быть bin или nc (netcdf). Если этот параметр отсутствует, формат определяется по расширению выходного файла (.bin или .nc).
+ velout - выходной файл компонентов скорости с прореженной сеткой. При отсутствии этот файл не пишется.
+ veloutformat - то же самое, что и outformat, только для velout файла.
+ shiftx, shifty - сдвиг начальной точки прореженной сетки относительно исходной. По умолчанию оба параметра равны нулю.
+ skipx, skipy - шаг по исходной сетке. Если равен 1, то выводится каждая точка исходной сетки, если 2, то каждая вторая и т. д. По умолчанию оба параметра равны 1.
+ stpout - выходной файл с особыми точками. При отсутствии этот файл не пишется.
+ stpoutformat - формат файла с особыми точками. Смотри outformat.
+ compress - число от 0 до 9, степень сжатия netcdf файлов. По умолчанию - 3.
+ obfuscate - если параметр присутствует, не записывать в выходные файлы метаинформацию.
+int - интерполирует данные для точек траекторий
+ var или vars - список разделённых запятыми физических величин. См. описание этого параметра для action tsc.
+ out - выходной файл.
+ input - входной файл с траекториями.
+genintfile - делает файл со скоростями для расчёта лагранжевых траекторий.
+ time, timeb, timee - смотри описание tsc.
+ lonb, lone, latb, late - смотри описание tsc.
+ geostrophic. Смотри описание этого параметра для action=uv.
+ out - выходной файл.
+grad - Работает только с источником TSCDATA. Берётся выходной файл от action=tsc и для всех переменных в нём рассчитывает градиенты по методу Белкина.
+ out - выходной файл.
+ compress - число от 0 до 9, степень сжатия netcdf файлов. По умолчанию - 3.
+ $var_min, $var_max, $var_log, где $var - имя переменной. Управляют первым этапом, нормализацией значений переменных. Значения переменных приводятся к интервалу 0-65534, значения меньше $var_min становятся нулём, больше $var_max - 65534. Если $var_min и/или $var_max не указаны, они выставляются в минимальное/максимальное значение для данной переменной. Если задана $var_log (с любым значением или вообще без него), то перед нормализацией берётся логарифм от значения переменной, в этом случае $var_min и $var_max должны быть больше нуля.
+mirror - зеркалирование данных в локальное хранилище. Своих параметров нет, они определяются источником.
+
+Источники
+AVISO - спутниковые данные для уровня моря и геострофических скоростей. Для этого источника доступно выделение отдельного региона.
+ dataset - DT, NRT, EckmanDT или EckmanNRT. По умолчанию DT.
+ layer или depth. Для EckmanDT и EckmanNRT доступны два слоя. Подробнее смотри описание этих параметров для источника NEMO.
+AVISOLOCAL - спутниковые данные для уровня моря и геострофических скоростей, локальное зеркало. Для этого источника доступно выделение отдельного региона. Других параметров нет.
+BINFILE - файл со скоростями для расчёта лагранжевых траекторий. Регион не нужен, всегда используется вся область файла.
+ dataset - путь или идентификатор файла.
+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.
+NEMO - данные реанализа по модели NEMO. Нужно указывать регион.
+ dataset - DT, NRT, NRT6 (глобальные), BALTICDT, BALTICNRT, BALTICNRT1, BLKSEADT, BLKSEANRT, MEDSEADT, MEDSEANRT, BISCDT, BISCNRT, ENWSDT, ENWSNRT (региональные). По умолчанию DT.
+ layer или depth - выбор номера слоя или глубины в метрах. Если указаны оба параметра, приоритет имеет depth. Если указан depth, выбирается слой, ближайший к заданной глубине. Нумерация слоёв начинается с нуля и идёт сверху вниз. Если не указаны оба параметра, используется нулевой слой (поверхность). Эти параметры игнорируются для некоторых переменных, например, ssh или mld.
+NEMOBIO - модель распространения биогенов по данным реанализа NEMO. Нужно указывать регион.
+ dataset - DT или NRT. По умолчанию DT.
+ layer или depth - выбор слоя или глубины. Подробнее смотри описание этих параметров для источника NEMO.
+VYLET - данные лагранжева моделирования с помощью vylet. Регион берётся из файла с данными.
+ dataset - имя файла с данными.
+TSCDATA - данные от tsc. Используется действием grad.
+ dataset - имя выходного файла от tsc в netcdf формате.
+GRIDFILE - данные для скорости в виде двух файлов в Surfer 6 grd формате. Время не используется и всегда равно началу эпохи.
+ u - файл с u-компонентой
+ v - файл с v-компонентой
+ unit - единица измерения. По умолчанию - cm/s
+COPERNICUS - источник для зеркалирования данных с серверов copernicus.
+ product - название продукта
+ dataset - названия наборов данных в конкретном продукте, через запятую. Если отсутствует, будут зеркалироваться все наборы данных продукта.
+ filter - регулярное выражение, определяющее имена файлов для зеркалирования. По умолчанию - ".*". Только для действия mirror.
+ Следующие параметры влияют на выбор подмножества данных и работают для действий tsc и info.
+ lonb, lone, latb, late - регион, для которого берутся данные. В случае отсутствия используются минимумы/максимумы из данных.
+ layer, layers, depth, depths - выбор горизонта/горизонтов. Если отсутствуют все эти параметры, используется слой с нулевым номером (приповерхностный). Если присутствуют два или более параметра, то использует параметр с большим приоритетом. layer > depth > layers > depths.
+  layer - номер горизонта, ноль соответствует поверхности,
+  layers - номер горизонта или диапазон горизонтов через двоеточие или слово ALL (all) для выбора всех горизонтов,
+  depth - глубина, выбирается горизонт, ближайший к этой глубине,
+  depths - глубина или диапазон глубин, разделённых двоеточием, горизонты выбираются так же, как и для параметра depth.
+ time, timeb, timee, timefilt - параметры, управляющие выбором моментов времени. Использовать можно либо time, либо тройку timeb, timee, timefilt.
+  time - фиксированое время, регулярное выражение или одно из значений BEGIN, BEG, FIRST (соответствуют минимальному времени, для которого доступны данные) или END, LAST (соответствуют максимальному времени, для которого доступны данные). В случае фиксированого времени, выбирается ближайшее доступное время. Если time регулярное выражение, то это синоним timefilt.
+  timeb, timee, timefilt - выбираются времена в интервале timeb : timee, удовлетворяющие регулярному выражению timefilt. Если timefilt отсутствует, используются все времена в интервале. По умолчанию timeb=BEG, timee=END.
+ var, vars - синонимы, var имеет больший приоритет. Список переменных, разделённых запятыми. Базовый список можно посмотреть в выводе действия info.

genusageh

@@ -0,0 +1,3 @@
+#!/bin/bash
+
+sed -e 's/"/\\"/g;s/$/\\n"/;s/^/"/' "$1" > "$2"

include/Adapter.h

@@ -0,0 +1,152 @@
+#pragma once
+#include "Data3D.h"
+#include "ParseArgs.h"
+#include "mdatetime.h"
+#include "merrors.h"
+#include "nczarr.h"
+#include <functional>
+
+using michlib::MDateTime;
+using michlib::RetVal;
+
+class Adapter
+{
+ public:
+ // Информация о блоке чтения, имена и индексы размерностей
+ struct ReadInfo
+ {
+  MString xdname, ydname, zdname, tdname;
+  size_t  xb, xe;
+  size_t  yb, ye;
+  size_t  zb, ze;
+
+  auto operator<=>(const ReadInfo&) const = default;
+ };
+
+ // Информация о переменной
+ struct VInfo
+ {
+  MString                          name;                         // Внутреннее имя в файле
+  MString                          unit, stname, lname, comment; // Единица измерения в файле, стандартное имя, человеческое имя, комментарий
+  MString                          targetunit;                   // Единица измерения прочитанных данных
+  std::shared_ptr<struct ReadInfo> rinfo;                        // Информация о блоке чтения
+  real                             scale;                        // Масштаб и
+  real                             offset;                       // смещение величины в файле
+  real                             fillval;                      // _FillVal
+ };
+
+ // Полная информация о переменной
+ struct VarInfo
+ {
+  struct VInfo info; // Информация о переменной в файле
+
+  std::shared_ptr<Projection> proj; // Горизонтальная структура прочитанных данных
+  std::shared_ptr<Vertical>   vert; // Вертикальная структура прочитанных данных
+
+  // Функции чтения двумерного и трёхмерного блоков
+  std::function<RetVal<std::shared_ptr<Data2D>>(const Adapter&, const void*, std::shared_ptr<Projection>, size_t)>                            Read2D; // = Def2DReader;
+  std::function<RetVal<std::shared_ptr<Data3D>>(const Adapter&, const void*, std::shared_ptr<Projection>, std::shared_ptr<Vertical>, size_t)> Read3D; // = Def3DReader;
+ };
+
+ // Запись в таблице файлов
+ struct TimeRow
+ {
+  size_t  beg, end; // Начальный и конечный временные индексы
+  MString file;     // Имя файла, в котором хранятся данные для этого интервала времени
+ };
+
+ enum class ZarrMethod
+ {
+  NC,
+  ZARR,
+  MNC,
+  MZARR
+ };
+
+ Adapter(): Adapter(ZarrMethod::NC) {};
+ Adapter(ZarrMethod m): method(m) {}
+ Adapter(Adapter&&)            = default;
+ Adapter& operator=(Adapter&&) = default;
+
+ template<class T>
+ requires std::same_as<std::vector<MDateTime>, std::decay_t<T>>
+ Error SetTimes(T&& t, const CLArgs& args, michlib_internal::ParameterListEx& pars)
+ {
+  times   = std::forward<T>(t);
+  curtime = times.size();
+  return FilterTimes(args, pars);
+ }
+
+ template<class T>
+ requires std::same_as<std::vector<TimeRow>, std::decay_t<T>>
+ void SetTimeTable(T&& t)
+ {
+  timetable = std::forward<T>(t);
+ }
+
+ void SetTitle(const MString& t) { title = t; }
+
+ void Add2DVariable(const MString& name, struct VInfo&& vinfo, std::shared_ptr<Projection> proj, decltype(VarInfo().Read2D) func = Def2DReader)
+ {
+  vars[name] = {.info = std::move(vinfo), .proj = proj, .Read2D = func, .Read3D = {}};
+ }
+
+ void Add3DVariable(const MString& name, struct VInfo&& vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, decltype(VarInfo().Read3D) func = Def3DReader)
+ {
+  vars[name] = {.info = std::move(vinfo), .proj = proj, .vert = vert, .Read2D = {}, .Read3D = func};
+ }
+
+ 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; }
+ const auto& Title() const { return title; }
+ const auto& Var(const MString& name) const { return vars.at(name); }
+
+ std::vector<MDateTime> IndexedTimes() const
+ {
+  std::vector<MDateTime> out(tindexes.size());
+  for(size_t i = 0; i < out.size(); i++) out[i] = times[tindexes[i]];
+  return out;
+ }
+
+ bool HasVar(const MString& name) const { return vars.contains(name); }
+
+ bool Is2D(const MString& name) const { return HasVar(name) && vars.at(name).Read2D; }
+ bool Is3D(const MString& name) const { return HasVar(name) && vars.at(name).Read3D; }
+
+ private:
+ // Функция чтения двумерных данных по умолчанию
+ static RetVal<std::shared_ptr<Data2D>> Def2DReader(const Adapter& ad, const void* vinfo, std::shared_ptr<Projection> proj, size_t it);
+ // Функция чтения трёхмерных данных по умолчанию
+ static RetVal<std::shared_ptr<Data3D>> Def3DReader(const Adapter& ad, const void* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it);
+
+ // Создание индексов tindexes сообразно параметрам командной строки
+ Error FilterTimes(const CLArgs& args, michlib_internal::ParameterListEx& pars);
+
+ MString                    title;     // Описание данных
+ std::unique_ptr<NCZarr>    nc;        // Текущий файл, с которого происходит чтение
+ MString                    curfile;   // Его имя
+ std::map<MString, VarInfo> vars;      // Переменные, которые можно читать
+ ZarrMethod                 method;    // Метод доступа к данным файла
+ std::vector<TimeRow>       timetable; // Таблица файлов, каждый файл содержит даннные для соответствующего временного интервала
+ size_t                     curtime;   // Номер текущей записи в таблице
+ size_t                     curit;     // Смещение текущего времени в текущем файле
+ std::vector<MDateTime>     times;     // Моменты времени, для которых доступны данные
+ std::vector<size_t>        tindexes;  // Индексы в times, для которых запрашиваются данные
+
+ // Кэши прочитаннх данных
+ std::unique_ptr<std::map<MString, std::shared_ptr<Data2D>>> cache2D;
+ std::unique_ptr<std::map<MString, std::shared_ptr<Data3D>>> cache3D;
+
+ // Проверка корректности запроса переменной и времени, подготовка адаптера к чтению данных
+ Error ReadCheck(const MString& var, size_t it);
+
+ public:
+ RetVal<std::shared_ptr<Data2D>> Read2D(const MString& var, size_t it);
+ RetVal<std::shared_ptr<Data3D>> Read3D(const MString& var, size_t it);
+};

include/CF.h

@@ -0,0 +1,80 @@
+#pragma once
+#include "Adapter.h"
+
+class CF
+{
+ public:
+ struct AxisInfo
+ {
+  MString tvarname, zvarname, yvarname, xvarname;
+  MString tdimname, zdimname, ydimname, xdimname;
+  uint    tdims = 0, zdims = 0, ydims = 0, xdims = 0;
+ };
+
+ // Get time variable name
+ static MString GetTVarName(const NCZarr& nc);
+
+ // Get info about dimensions
+ static struct AxisInfo GetAxisInfo(const NCZarr& nc);
+
+ // Parse reference date and time step from time description
+ static std::tuple<MDateTime, time_t, bool> ParseRefDate(const MString& refdate);
+
+ // Read time moments from file
+ static RetVal<std::vector<MDateTime>> ReadTimes(const NCZarr& nc, const MString& tname);
+
+ // Check, if vector data lays on uniform grid
+ static bool CheckUniform(std::ranges::random_access_range auto&& data, real tolerance = 1e-8)
+ requires std::convertible_to<std::ranges::range_value_t<decltype(data)>, real>
+ {
+  if(data.size() < 2) return false;
+  if(data.size() == 2) return true;
+
+  real step = data[1] - data[0];
+
+  for(size_t i = 1; i < data.size() - 1; i++)
+   if(michlib::Abs(1.0 - (data[i + 1] - data[i]) / step) > tolerance) return false;
+  return true;
+ }
+
+ // Check, if vector time data lays on uniform grid
+ static bool CheckUniform(std::ranges::random_access_range auto&& data)
+ requires std::convertible_to<std::ranges::range_value_t<decltype(data)>, MDateTime>
+ {
+  if(data.size() < 2) return false;
+  if(data.size() == 2) return true;
+
+  auto step = data[1] - data[0];
+
+  for(size_t i = 1; i < data.size() - 1; i++)
+   if(data[i + 1] - data[i] != step) return false;
+  return true;
+ }
+
+ static Error FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad, std::unique_ptr<NCZarr>&& pnc);
+
+ static MString StName2Name(const MString& stname)
+ {
+  if(stname == "sea_water_potential_temperature") return "ptemp";
+  if(stname == "sea_water_temperature") return "temp";
+  if(stname == "sea_water_salinity") return "sal";
+  if(stname == "ocean_mixed_layer_thickness_defined_by_sigma_theta") return "mld";
+  if(stname == "sea_surface_height_above_geoid") return "ssh";
+  if(stname == "sea_surface_elevation") return "ssh";
+  if(stname == "eastward_sea_water_velocity") return "u";
+  if(stname == "northward_sea_water_velocity") return "v";
+  if(stname == "upward_sea_water_velocity") return "w";
+  if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "ugs";
+  if(stname == "surface_geostrophic_northward_sea_water_velocity") return "vgs";
+  if(stname == "mass_concentration_of_chlorophyll_a_in_sea_water") return "chl";
+  if(stname == "mole_concentration_of_nitrate_in_sea_water") return "NO3";
+  if(stname == "net_primary_production_of_biomass_expressed_as_carbon_per_unit_volume_in_sea_water") return "prprod";
+  if(stname == "mole_concentration_of_phytoplankton_expressed_as_carbon_in_sea_water") return "Cchl";
+  if(stname == "mole_concentration_of_phosphate_in_sea_water") return "PO4";
+  if(stname == "mole_concentration_of_silicate_in_sea_water") return "Si";
+  if(stname == "mole_concentration_of_dissolved_molecular_oxygen_in_sea_water") return "O2";
+  return "";
+ }
+
+ //static RetVal<std::shared_ptr<Data2D>> PTemp3DReader(const Adapter& ad, const void* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it);
+};

include/Data3D.h

@@ -0,0 +1,740 @@
+#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(); }
+};

include/actiondep.h

@@ -1,6 +1,7 @@
 #pragma once
 #include "ParseArgs.h"
 #include "basedata.h"
+#include "conscolors.h"
 #include "mdatetime.h"
 #include "mregex.h"
 #include "traits.h"

include/conscolors.h

@@ -0,0 +1,67 @@
+#pragma once
+#include "merrors.h"
+
+using michlib::int_cast;
+using michlib::MString;
+using michlib::uint2;
+
+class C
+{
+ protected:
+ static bool tty;
+ MString     command;
+
+ template<class... T> static MString Create(uint2 c, T... args) { return MString(c) + ";" + Create(args...); }
+
+ static MString Create(uint2 c) { return MString(c); }
+
+ template<class... T> C(T... args): command("\033[" + Create(args...) + "m") {}
+
+ public:
+ static void Init() { tty = isatty(michlib_internal::fd_std) == 1; }
+ static void Init(bool c) { tty = c; }
+
+ enum
+ {
+  BLACK = 30,
+  RED,
+  GREEN,
+  YELLOW,
+  BLUE,
+  MAGENTA,
+  CYAN,
+  WHITE,
+  BBLACK = 90,
+  BRED,
+  BGREEN,
+  BYELLOW,
+  BBLUE,
+  BMAGENTA,
+  BCYAN,
+  BWHITE
+ };
+
+ C(): C(0) {}
+
+ operator MString()
+ {
+  if(tty)
+   return command;
+  else
+   return "";
+ }
+};
+
+class B: public C
+{
+ public:
+ B(): C(49) {}
+ B(decltype(C::RED) c): C(int_cast<uint2>(c) + 10) {}
+};
+
+class F: public C
+{
+ public:
+ F(): C(39) {}
+ F(decltype(C::RED) c): C(int_cast<uint2>(c)) {}
+};

include/copcat.h

@@ -36,9 +36,15 @@ class CopernicusCatalog
  // URL of product
  RetVal<MString> ProductURL(const MString& prod) const;
 
+ // Title of product
+ RetVal<MString> ProductTittle(const MString& prod) const;
+
  // URL of dataset
  RetVal<MString> DatasetURL(const MString& prod, const MString& dataset) const;
 
+ // Title of dataset
+ RetVal<MString> DatasetTitle(const MString& prod, const MString& dataset) const;
+
  // URL of native data (files) in dataset
  RetVal<MString> DatasetNativeURL(const MString& prod, const MString& dataset) const { return AssetURL(prod, dataset, "native"); }
 

include/ncfilew.h

@@ -1,4 +1,5 @@
 #pragma once
+#include "Adapter.h"
 #include "MString.h"
 #include "actiondep.h"
 #include "traits.h"
@@ -386,6 +387,14 @@ class NCFileW: public NCFileWBase
  public:
  static constexpr auto Fill() { return fill; }
 
+ michlib::Error AddTimeData(const std::vector<MDateTime>& times);
+
+ michlib::Error Create(const MString& name, const Adapter& ad, const michlib_internal::ParameterListEx& pars, int compression);
+
+ template<class D>
+ 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)
  {
   if(type != UNKNOWN) return "File already created";

include/traits.h

@@ -1,4 +1,5 @@
 #pragma once
+#include "Adapter.h"
 #include "ParseArgs.h"
 #include <concepts>
 
@@ -6,32 +7,29 @@ class BaseParameters;
 
 using michlib::real;
 
+template<class T>
+concept AdapterSupported = requires(T t, const CLArgs& args, michlib_internal::ParameterListEx& pars) {
+ { t.GetAdapter(args, pars) } -> std::convertible_to<RetVal<Adapter>>;
+};
+
 template<class T>
 concept InfoSupported = requires(T t) {
- {
-  t.Info()
- } -> std::convertible_to<MString>;
+ { t.Info() } -> std::convertible_to<MString>;
 };
 
 template<class T>
 concept HasDefVars = requires(T t) {
- {
-  t.DefaultVars()
- } -> std::convertible_to<MString>;
+ { t.DefaultVars() } -> std::convertible_to<MString>;
 };
 
 template<class T>
 concept ParametersRequiredRegion = requires(T t, michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& r) {
- {
-  t.Parameters(pars, args, r).first
- } -> std::convertible_to<const BaseParameters*>;
+ { t.Parameters(pars, args, r).first } -> std::convertible_to<const BaseParameters*>;
 };
 
 template<class T>
 concept ParametersNotRequiredRegion = requires(T t, michlib_internal::ParameterListEx& pars, const CLArgs& args) {
- {
-  t.Parameters(pars, args).first
- } -> std::convertible_to<const BaseParameters*>;
+ { t.Parameters(pars, args).first } -> std::convertible_to<const BaseParameters*>;
 };
 
 template<class T>
@@ -39,40 +37,26 @@ concept ParametersSupported = ParametersRequiredRegion<T> || ParametersNotRequir
 
 template<class T>
 concept ReadPSupported = requires(T t, const MString& vname, std::map<MString, typename T::Data>& cache, const BaseParameters* ip, size_t i) {
- {
-  t.Read(vname, cache, ip, i)
- } -> std::same_as<bool>;
+ { t.Read(vname, cache, ip, i) } -> std::same_as<bool>;
 };
 template<class T>
 concept ReadSupported = requires(T t, const MString& vname, std::map<MString, typename T::Data>& cache, size_t i) {
- {
-  t.Read(vname, cache, i)
- } -> std::same_as<bool>;
+ { t.Read(vname, cache, i) } -> std::same_as<bool>;
 };
 
 template<class T>
 concept IsUVData = requires(T t) {
- {
-  t.U(0)
- } -> std::convertible_to<real>;
- {
-  t.V(0)
- } -> std::convertible_to<real>;
- {
-  t.U2(0)
- } -> std::convertible_to<real>;
+ { t.U(0) } -> std::convertible_to<real>;
+ { t.V(0) } -> std::convertible_to<real>;
+ { t.U2(0) } -> std::convertible_to<real>;
 };
 
 namespace internal
 {
 template<class T>
 concept HaveXYStep = requires(T t) {
- {
-  t.XStep()
- } -> std::convertible_to<real>;
- {
-  t.YStep()
- } -> std::convertible_to<real>;
+ { t.XStep() } -> std::convertible_to<real>;
+ { t.YStep() } -> std::convertible_to<real>;
 };
 
 template<typename...> using void_ = void;
@@ -91,9 +75,7 @@ template<class D> struct GetReadType_<D, false>
 
 template<class T>
 concept HaveDisable = requires(T t) {
- {
-  T::disabledactions
- } -> std::convertible_to<const char*>;
+ { T::disabledactions } -> std::convertible_to<const char*>;
 };
 
 consteval bool cmpspace(const char* s1, const char* s2)
@@ -130,70 +112,42 @@ template<class D> using ReadType = internal::GetReadType_<D, ReadPSupported<D> |
 
 template<class T>
 concept ReadIsGrid = requires {
- {
-  std::declval<ReadType<T>>().XStep()
- } -> std::convertible_to<real>;
- {
-  std::declval<ReadType<T>>().YStep()
- } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().XStep() } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().YStep() } -> std::convertible_to<real>;
 };
 
 template<class T>
 concept ReadIs2DGeoRectArray = requires {
- {
-  std::declval<ReadType<T>>().Ix2Lon(0)
- } -> std::convertible_to<real>;
- {
-  std::declval<ReadType<T>>().Iy2Lat(0)
- } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Ix2Lon(0) } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Iy2Lat(0) } -> std::convertible_to<real>;
 };
 
 template<class T>
 concept ReadIs2DGeoArray = requires {
- {
-  std::declval<ReadType<T>>().Lon(0, 0)
- } -> std::convertible_to<real>;
- {
-  std::declval<ReadType<T>>().Lat(0, 0)
- } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Lon(0, 0) } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Lat(0, 0) } -> std::convertible_to<real>;
 };
 
 template<class T>
 concept ReadIs1DGeoArray = requires {
- {
-  std::declval<ReadType<T>>().Lon(0)
- } -> std::convertible_to<real>;
- {
-  std::declval<ReadType<T>>().Lat(0)
- } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Lon(0) } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Lat(0) } -> std::convertible_to<real>;
 };
 
 template<class T>
 concept ReadIs2DXYRectArray = requires {
- {
-  std::declval<ReadType<T>>().Ix2X(0)
- } -> std::convertible_to<real>;
- {
-  std::declval<ReadType<T>>().Iy2Y(0)
- } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Ix2X(0) } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Iy2Y(0) } -> std::convertible_to<real>;
 };
 
 template<class T>
 concept ReadIs2DXYArray = requires {
- {
-  std::declval<ReadType<T>>().X(0, 0)
- } -> std::convertible_to<real>;
- {
-  std::declval<ReadType<T>>().Y(0, 0)
- } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().X(0, 0) } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Y(0, 0) } -> std::convertible_to<real>;
 };
 
 template<class T>
 concept ReadIs1DArray = requires {
- {
-  std::declval<ReadType<T>>().X(0)
- } -> std::convertible_to<real>;
- {
-  std::declval<ReadType<T>>().Y(0)
- } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().X(0) } -> std::convertible_to<real>;
+ { std::declval<ReadType<T>>().Y(0) } -> std::convertible_to<real>;
 };

include/udunitscpp.h

@@ -0,0 +1,94 @@
+#pragma once
+#include "MString.h"
+#include <udunits2.h>
+
+using michlib::MString;
+
+class UtUnit
+{
+ static ut_system* system;
+ static size_t     count;
+
+ ut_unit* unit;
+
+ class ConverterRAII
+ {
+  cv_converter* cnv;
+
+  ConverterRAII()                     = delete;
+  ConverterRAII(const ConverterRAII&) = delete;
+
+  public:
+  ConverterRAII(ConverterRAII&&) = default;
+
+  ConverterRAII(const ut_unit* fr, const ut_unit* to): cnv(ut_get_converter(fr, to))
+  {
+   /*
+   auto* frc = ut_clone(fr);
+   auto* toc = ut_clone(to);
+   cnv       = ut_get_converter(frc, toc);
+   if(frc != nullptr) ut_free(frc);
+   if(toc != nullptr) ut_free(toc);
+   */
+  }
+  ~ConverterRAII()
+  {
+   if(cnv != nullptr) cv_free(cnv);
+  }
+
+  explicit operator bool() const { return cnv != nullptr; }
+  float    operator()(float v) const { return cnv == nullptr ? v : cv_convert_float(cnv, v); }
+  double   operator()(double v) const { return cnv == nullptr ? v : cv_convert_float(cnv, v); }
+ };
+
+ static void ReadSystem()
+ {
+  if(system == nullptr)
+  {
+   auto oldhandler = ut_set_error_message_handler(ut_ignore);
+   system          = ut_read_xml(nullptr);
+   ut_set_error_message_handler(oldhandler);
+  }
+ }
+
+ static void DestroySystem()
+ {
+  if(system != nullptr) ut_free_system(system);
+  system = nullptr;
+ }
+
+ UtUnit()              = delete;
+ UtUnit(const UtUnit&) = delete;
+
+ public:
+ UtUnit(UtUnit&&) = default;
+
+ UtUnit(const MString& desc)
+ {
+  ReadSystem();
+  unit = ut_parse(system, desc.Buf(), UT_UTF8);
+  if(unit != nullptr) count++;
+ }
+
+ ~UtUnit()
+ {
+  if(unit != nullptr)
+  {
+   count--;
+   ut_free(unit);
+   unit = nullptr;
+  }
+  if(count == 0) DestroySystem();
+ }
+
+ explicit operator bool() const { return unit != nullptr; }
+
+ bool operator==(const UtUnit& u) const { return ut_compare(unit, u.unit) == 0; }
+
+ bool Compatible(const UtUnit& u) const { return ut_are_convertible(unit, u.unit) != 0; }
+
+ auto Converter(const UtUnit& u) const { return ConverterRAII{unit, u.unit}; }
+};
+
+inline bool Compatible(const UtUnit& u1, const UtUnit& u2) { return u1.Compatible(u2); }
+inline auto Converter(const UtUnit& fr, const UtUnit& to) { return fr.Converter(to); }

michlib

@@ -1 +1 @@
-Subproject commit f380988909fadd7f42c7ce09288c4ff3198ca318
+Subproject commit 1e53f661ba1697a11c74205452eae4e598112ae9

sources/COPERNICUS.cpp

@@ -1,8 +1,11 @@
 #define MICHLIB_NOSOURCE
 #include "COPERNICUS.h"
+#include "CF.h"
 #include "mirrorfuncs.h"
 #include <libxml/parser.h>
 #include <libxml/tree.h>
+#include <numeric>
+#include <valarray>
 
 using michlib::GPL;
 
@@ -128,7 +131,7 @@ Error COPERNICUSData::Mirror(const CLArgs& args) const
 
  std::vector<MString> dsets;
  if(args.contains("dataset"))
-  dsets.push_back(args.at("dataset"));
+  dsets = args.at("dataset").Split(":,");
  else
  {
   auto dlist = cat.DatasetList(prod);
@@ -232,3 +235,123 @@ Error COPERNICUSData::Mirror(const CLArgs& args) const
 
  return Error();
 }
+
+RetVal<Adapter> COPERNICUSData::GetAdapter(const CLArgs& args, michlib_internal::ParameterListEx& pars) const
+{
+ static const MString pref = "COPERNICUSData::GetAdapter";
+
+ bool debug = args.contains("debug");
+
+ Adapter out(Adapter::ZarrMethod::MZARR);
+
+ GPL.UsePrefix("NEMO");
+
+ // Get link for dataset
+ MString              product;
+ std::vector<MString> dset;
+ MString              dataset;
+ {
+  if(args.contains("product")) // Get product and dataset from command line
+  {
+   product = args.at("product");
+   if(args.contains("dataset")) dset = args.at("dataset").Split(":,");
+   dataset = product + " " + args.at("dataset");
+  }
+  else // Get product and dataset from command line
+  {
+   // Try one url
+   dataset     = args.contains("dataset") ? args.at("dataset") : "DT";
+   MString url = GPL.ParameterSValue(dataset + "_URL", "");
+   if(url.Exist())
+   {
+    auto list = url.Split(":");
+    product   = list[0];
+    for(size_t i = 1; i < list.size(); i++) dset.push_back(list[i]);
+   }
+   else // Multiple urls
+   {
+    size_t i = 1;
+    while(true)
+    {
+     MString url = GPL.ParameterSValue(dataset + "_URL" + i, "");
+     if(url.Exist())
+     {
+      // Split url on product and dataset
+      auto words = url.Split(":");
+      if(words.size() != 2) return {pref, "Invalid url for dataset " + dataset + ": " + url};
+      if(product.Exist() && product != words[0]) return {pref, "Product for dataset " + dataset + " was already defined: " + url};
+      product = words[0];
+      dset.push_back(words[1]);
+     }
+     else
+      break;
+     i++;
+    }
+   }
+   if(!product.Exist()) return {pref, "Can't find product for the dataset " + dataset};
+   pars.AddParameter("dataset", dataset);
+  }
+
+  if(dset.size() == 0) // Get datasets from catalog
+  {
+   CopernicusCatalog cat;
+   auto              dlist = cat.DatasetList(product);
+   if(!dlist) return dlist.Add(pref, "Can't get list of datasets for product " + product);
+   dset = dlist.Value();
+  }
+ }
+
+ pars.AddParameter("product", product);
+ pars.AddParameter("datasets", std::accumulate(dset.begin() + 1, dset.end(), dset.front(), [](const MString& a, const MString& b) { return a + "," + b; }));
+
+ // Open data
+ std::unique_ptr<NCZarr> pnc(new NCZarr);
+ {
+  auto ret = pnc->OpenMultiZarr(product, dset);
+  if(!ret) return ret.Add(pref, "Can't open dataset " + dataset);
+ }
+
+ // Title
+ {
+  MString           title;
+  CopernicusCatalog cat;
+  auto              ret = cat.DatasetTitle(product, dset[0]);
+  if(!ret) return ret.Add(pref, "Can't get title for dataset " + dataset);
+  title = ret.Value();
+  pars.AddParameter("title", title);
+  out.SetTitle(title);
+ }
+
+ // Filename for adapter
+ MString fname = product + ":" + dset[0];
+ for(size_t i = 1; i < dset.size(); i++) fname += "," + dset[i];
+
+ 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))};
+ if(debug) pars.AddParameter("File name", fname);
+
+ // Read and set times
+ {
+  auto ret = CF::ReadTimes(*pnc, tvar);
+  if(!ret) return ret.Add(pref, "Can't read time values from the dataset " + dataset);
+  auto ret1 = out.SetTimes(ret.Value(), args, pars);
+  if(!ret1) return ret1.Add(pref, "Can't filter time values from the dataset " + dataset);
+ }
+
+ // Timetable
+ {
+  std::vector<Adapter::TimeRow> table{1};
+  table[0].beg  = 0;
+  table[0].end  = out.Times().size() - 1;
+  table[0].file = fname;
+  out.SetTimeTable(std::move(table));
+ }
+
+ {
+  auto ret = CF::FillAdapterFromCF(args, pars, out, std::move(pnc));
+  if(!ret) return ret.Add(pref, "Can't get adapter");
+ }
+
+ return out;
+}

sources/COPERNICUS.h

@@ -1,10 +1,13 @@
 #pragma once
+#include "Adapter.h"
 #include "ParseArgs.h"
 #include "copcat.h"
 #include "mdatetime.h"
 
+using michlib::DetGeoDomain;
 using michlib::MDateTime;
 using michlib::MString;
+using michlib::ToGeoDomain;
 
 class COPERNICUSData
 {
@@ -18,4 +21,7 @@ class COPERNICUSData
 
  // Main mirror function
  Error Mirror(const CLArgs& args) const;
+
+ // Adapter
+ RetVal<Adapter> GetAdapter(const CLArgs& args, michlib_internal::ParameterListEx& pars) const;
 };

src/Adapter.cpp

@@ -0,0 +1,363 @@
+#define MICHLIB_NOSOURCE
+#include "Adapter.h"
+
+Error Adapter::FilterTimes(const CLArgs& args, michlib_internal::ParameterListEx& pars)
+{
+ static const MString pref = "Adapter::FilterTimes";
+
+ if(args.contains("time") && (args.contains("timeb") || args.contains("timee"))) return {pref, "Time must be set via time parameter or timeb and timee parameter but not via both"};
+ if(args.contains("time") && args.contains("timefilt")) return {pref, "Time filter must be set via time parameter or timefilter parameter but not via both"};
+
+ MDateTime b, e;
+
+ MString tb = args.contains("timeb") ? args.at("timeb") : "BEG";
+ MString te = args.contains("timee") ? args.at("timee") : "END";
+ MString tf = args.contains("timefilt") ? args.at("timefilt") : "";
+
+ if(args.contains("time"))
+ {
+  MDateTime time;
+  MString   t = args.at("time");
+
+  if(t == "BEGIN" || t == "BEG" || t == "FIRST")
+   t = times.front().ToTString();
+  else if(t == "END" || t == "LAST")
+   t = times.back().ToTString();
+  else if(time.FromString(t)) // Time, not regex
+  {
+   if(time <= times.front())
+    t = times.front().ToTString();
+   else if(time >= times.back())
+    t = times.back().ToTString();
+   else
+    for(size_t i = 0; i < times.size() - 1; i++)
+    {
+     if(time >= times[i] && time <= times[i + 1])
+     {
+      t = times[(time - times[i] <= times[i + 1] - time) ? i : (i + 1)].ToTString();
+      break;
+     }
+    }
+   tb = te = t;
+  }
+  else // Regular expression
+   tf = t;
+ }
+
+ if(tb == "BEGIN" || tb == "BEG" || tb == "FIRST")
+  b = times.front();
+ else if(tb == "END" || tb == "LAST")
+  b = times.back();
+ else
+  b.FromString(tb);
+
+ if(te == "BEGIN" || te == "BEG" || te == "FIRST")
+  e = times.front();
+ else if(te == "END" || te == "LAST")
+  e = times.back();
+ else
+  e.FromString(te);
+
+ if(e < b) return {pref, "Timeb must not be greater then timee"};
+
+ michlib::RegExpSimple reg(tf.Buf());
+ bool                  filt = tf.Exist();
+
+ if(filt && reg.Compile() != 0) return {pref, "Bad regular expression: " + tf};
+
+ tindexes.clear();
+
+ for(size_t it = 0; it < times.size(); it++)
+  if(times[it] >= b && times[it] <= e && (!filt || reg.Match(times[it].ToTString().Buf()))) tindexes.push_back(it);
+
+ if(tindexes.size() == 0) return {pref, "No times with given criteria are found"};
+
+ if(tindexes.size() == 1)
+  pars.AddParameter("time", times[tindexes[0]].ToTString());
+ else
+ {
+  pars.AddParameter("timebegin", times[tindexes.front()].ToTString());
+  pars.AddParameter("timeend", times[tindexes.back()].ToTString());
+  if(filt) pars.AddParameter("timefilter", tf);
+ }
+
+ return {};
+}
+
+RetVal<std::shared_ptr<Data2D>> Adapter::Def2DReader(const Adapter& ad, const void* vinfo, std::shared_ptr<Projection> proj, size_t it)
+{
+ static const MString pref = "Adapter::Def2DReader";
+
+ auto v = michlib::pointer_cast<const struct VInfo*>(vinfo);
+
+ const struct ReadInfo* ri = v->rinfo.get();
+
+ const bool layered = ri->zdname.Exist() && ad.nc->HasDim(v->name, ri->zdname);
+
+ std::shared_ptr<Data2D> out(new Data2D(proj, {.unit = v->targetunit, .stname = v->stname, .lname = v->lname, .comment = v->comment, .fillval = v->fillval}));
+
+ auto trans = [scale = v->scale, offset = v->offset](auto raw) -> real { return raw * scale + offset; };
+ auto data  = [pdata = out.get()](size_t ix, size_t iy) -> real& { return (*pdata)(ix, iy); };
+
+ const UtUnit fr(v->unit), to(v->targetunit);
+ const auto   cnv      = Converter(fr, to);
+ auto         cnvtrans = [scale = v->scale, offset = v->offset, &cnv = std::as_const(cnv)](auto raw) -> real { return cnv(raw * scale + offset); };
+
+ const bool needsconvert = (v->unit.Exist() && v->targetunit.Exist() && fr && to && fr != to);
+ if(needsconvert && !cnv) return {pref, "Can't convert units from \"" + v->unit + "\" to \"" + v->targetunit + "\""};
+
+ MString xreq, yreq, zreq, treq;
+
+ yreq = ri->ydname + ":" + ri->yb + ":" + (ri->ye - ri->yb + 1);
+ zreq = layered ? (ri->zdname + ":" + ri->zb + ":1") : "";
+ treq = ri->tdname + ":" + it + ":1";
+
+ auto req = [&x = xreq, &y = yreq, &z = zreq, &t = treq]() -> MString { return x + ";" + y + ";" + (z.Exist() ? (z + ";") : "") + t; };
+
+ if(ri->xb < ri->xe)
+ {
+  xreq = ri->xdname + ":" + ri->xb + ":" + (ri->xe - ri->xb + 1);
+
+  auto ret = needsconvert ? ad.nc->Read(v->name, data, cnvtrans, req()) : ad.nc->Read(v->name, data, trans, req());
+  if(!ret) return ret.Add(pref, "Can't read variable " + v->name);
+ }
+ else
+ {
+  size_t nx = ad.nc->DimSize(v->name, ri->xdname);
+  {
+   xreq = ri->xdname + ":" + ri->xb + ":" + (nx - ri->xb + 1);
+
+   auto ret = needsconvert ? ad.nc->Read(v->name, data, cnvtrans, req()) : ad.nc->Read(v->name, data, trans, req());
+   if(!ret) return ret.Add(pref, "Can't read variable " + v->name);
+  }
+  {
+   size_t shift       = nx - ri->xb + 1;
+   auto   shifteddata = [pdata = out.get(), shift](size_t ix, size_t iy) -> real& { return (*pdata)(ix + shift, iy); };
+   xreq               = ri->xdname + ":0:" + (ri->xe - 1);
+
+   auto ret = needsconvert ? ad.nc->Read(v->name, shifteddata, cnvtrans, req()) : ad.nc->Read(v->name, shifteddata, trans, req());
+   if(!ret) return ret.Add(pref, "Can't read variable " + v->name);
+  }
+ }
+ return out;
+}
+
+RetVal<std::shared_ptr<Data3D>> Adapter::Def3DReader(const Adapter& ad, const void* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it)
+{
+ static const MString pref = "Adapter::Def3DReader";
+
+ auto v = michlib::pointer_cast<const struct VInfo*>(vinfo);
+
+ const struct ReadInfo* ri = v->rinfo.get();
+
+ const bool layered = ri->zdname.Exist() && ad.nc->HasDim(v->name, ri->zdname);
+ const bool invertz = ri->zb > ri->ze;
+ const bool onez    = ri->zb == ri->ze;
+
+ std::shared_ptr<Data3D> out(new Data3D(proj, vert, {.unit = v->targetunit, .stname = v->stname, .lname = v->lname, .comment = v->comment, .fillval = v->fillval}));
+
+ auto trans  = [scale = v->scale, offset = v->offset](auto raw) -> real { return raw * scale + offset; };
+ auto data   = [pdata = out.get(), invertz](size_t ix, size_t iy, size_t iz) -> real& { return (*pdata)(ix, iy, invertz ? (pdata->Nz() - iz - 1) : iz); };
+ auto data2D = [pdata = out.get()](size_t ix, size_t iy) -> real& { return (*pdata)(ix, iy, 0); };
+
+ const UtUnit fr(v->unit), to(v->targetunit);
+ const auto   cnv      = Converter(fr, to);
+ auto         cnvtrans = [scale = v->scale, offset = v->offset, &cnv = std::as_const(cnv)](auto raw) -> real { return cnv(raw * scale + offset); };
+
+ const bool needsconvert = (v->unit.Exist() && v->targetunit.Exist() && fr && to && fr != to);
+ if(needsconvert && !cnv) return {pref, "Can't convert units from \"" + v->unit + "\" to \"" + v->targetunit + "\""};
+
+ MString xreq, yreq, zreq, treq;
+
+ yreq = ri->ydname + ":" + ri->yb + ":" + (ri->ye - ri->yb + 1);
+ zreq = layered ? (ri->zdname + ":" + std::min(ri->zb, ri->ze) + ":" + (std::max(ri->zb, ri->ze) - std::min(ri->zb, ri->ze) + 1)) : "";
+ treq = ri->tdname + ":" + it + ":1";
+
+ auto req = [&x = xreq, &y = yreq, &z = zreq, &t = treq]() -> MString { return x + ";" + y + ";" + (z.Exist() ? (z + ";") : "") + t; };
+
+ if(ri->xb < ri->xe)
+ {
+  xreq = ri->xdname + ":" + ri->xb + ":" + (ri->xe - ri->xb + 1);
+
+  auto ret = needsconvert ? (onez ? ad.nc->Read(v->name, data2D, cnvtrans, req()) : ad.nc->Read(v->name, data, cnvtrans, req()))
+                          : (onez ? ad.nc->Read(v->name, data2D, trans, req()) : ad.nc->Read(v->name, data, trans, req()));
+  if(!ret) return ret.Add(pref, "Can't read variable " + v->name);
+ }
+ else
+ {
+  size_t nx = ad.nc->DimSize(v->name, ri->xdname);
+  {
+   xreq = ri->xdname + ":" + ri->xb + ":" + (nx - ri->xb + 1);
+
+   auto ret = needsconvert ? ad.nc->Read(v->name, data, cnvtrans, req()) : ad.nc->Read(v->name, data, trans, req());
+   if(!ret) return ret.Add(pref, "Can't read variable " + v->name);
+  }
+  {
+   size_t shift         = nx - ri->xb + 1;
+   auto   shifteddata   = [pdata = out.get(), shift, invertz](size_t ix, size_t iy, size_t iz) -> real& { return (*pdata)(ix + shift, iy, invertz ? (pdata->Nz() - iz - 1) : iz); };
+   auto   shifteddata2D = [pdata = out.get(), shift](size_t ix, size_t iy) -> real& { return (*pdata)(ix + shift, iy, 0); };
+   xreq                 = ri->xdname + ":0:" + (ri->xe - 1);
+
+   auto ret = needsconvert ? (onez ? ad.nc->Read(v->name, shifteddata2D, cnvtrans, req()) : ad.nc->Read(v->name, shifteddata, cnvtrans, req()))
+                           : (onez ? ad.nc->Read(v->name, shifteddata2D, trans, req()) : ad.nc->Read(v->name, shifteddata, trans, req()));
+   if(!ret) return ret.Add(pref, "Can't read variable " + v->name);
+  }
+ }
+ return out;
+}
+
+Error Adapter::ReadCheck(const MString& var, size_t it)
+{
+ static const MString pref = "Adapter::ReadCheck";
+
+ // Check if time is correct
+ if(it > times.size()) return {pref, MString("Time index ") + it + " is greater then number of times " + times.size()};
+
+ // Check presence of variable
+ if(!vars.contains(var)) return {pref, "Variable " + var + " not found"};
+
+ // Find and check entry in timetable
+ size_t timerow = timetable.size();
+ for(size_t i = 0; i < timetable.size(); i++)
+  if(it >= timetable[i].beg && it <= timetable[i].end)
+  {
+   timerow = i;
+   break;
+  }
+
+ if(timerow == timetable.size()) return {pref, "Time table does'nt contains row for time " + times[it].ToTString() + "/" + it};
+
+ // Time index inside file
+ curit = it - timetable[timerow].beg;
+
+ // Reopen data file if needed
+ if(timetable[timerow].file != curfile)
+ {
+  nc.reset(new NCZarr);
+  curfile = timetable[timerow].file;
+  switch(method)
+  {
+  case(ZarrMethod::NC):
+  {
+   auto ret = nc->OpenNC(curfile);
+   if(!ret) return ret.Add(pref, "OpenNC failed for file " + curfile);
+   break;
+  }
+  case(ZarrMethod::ZARR):
+  {
+   auto args = curfile.Split(":");
+   if(args.size() < 2 || args.size() > 3) return {pref, "Can't open Zarr file " + curfile + ": incorrect url"};
+   auto ret = (args.size() == 2) ? nc->OpenZarr(args[0], args[1]) : nc->OpenZarr(args[0], args[1], args[2].ToBool());
+   if(!ret) return ret.Add(pref, "OpenZarr failed for file " + curfile);
+   break;
+  }
+  case(ZarrMethod::MNC):
+  {
+   auto ret = nc->OpenMultiNC(curfile.Split(":"));
+   if(!ret) return ret.Add(pref, "OpenMultiNC failed for file " + curfile);
+   break;
+  }
+  case(ZarrMethod::MZARR):
+  {
+   auto args = curfile.Split(":");
+   if(args.size() < 2 || args.size() > 3) return {pref, "Can't open Zarr file " + curfile + ": incorrect url"};
+   auto ret = (args.size() == 2) ? nc->OpenMultiZarr(args[0], args[1].Split(",")) : nc->OpenMultiZarr(args[0], args[1].Split(","), args[2].ToBool());
+   if(!ret) return ret.Add(pref, "OpenMultiZarr failed for file " + curfile);
+   break;
+  }
+  }
+ }
+
+ {
+  const auto& v    = vars[var];
+  const auto* ri   = v.info.rinfo.get();
+  const auto& name = v.info.name;
+
+  const bool layered = ri->zdname.Exist() && nc->HasDim(name, ri->zdname);
+
+  if(layered)
+  {
+   //if(ri->ze < ri->zb) return {pref, "Internal error: ze < zb for variable " + var + "(" + name + ")" + " (" + ri->ze + " < " + ri->zb + ")"};
+   const size_t nz = nc->DimSize(name, ri->zdname);
+   if(ri->zb > nz - 1) return {pref, "Internal error: zb > nz - 1 for variable " + var + "(" + name + ")" + " (" + ri->zb + " > " + nz + " - 1)"};
+   if(ri->ze > nz - 1) return {pref, "Internal error: ze > nz - 1 for variable " + var + "(" + name + ")" + " (" + ri->ze + " > " + nz + " - 1)"};
+   if(v.vert && (ri->zb < ri->ze ? ri->ze - ri->zb : ri->zb - ri->ze) + 1 != v.vert->Nz())
+    return {pref, "Internal error: number of requested layers does'nt correspond parameters of vertical column for variable " + var + "(" + name + ")" + " (zb = " + ri->zb +
+                      ", ze = " + ri->ze + ", Nz = " + v.vert->Nz()};
+  }
+
+  if(ri->ye < ri->yb) return {pref, "Internal error: ye < yb for variable " + var + "(" + name + ")" + " (" + ri->ye + " < " + ri->yb + ")"};
+  const size_t nx = nc->DimSize(name, ri->xdname);
+  const size_t ny = nc->DimSize(name, ri->ydname);
+
+  if(ri->ye > ny - 1) return {pref, "Internal error: ye > ny - 1 for variable " + var + "(" + name + ")" + " (" + ri->ye + " > " + ny + " - 1)"};
+  if(ri->ye - ri->yb + 1 != v.proj->Ny())
+   return {pref, "Internal error: number of requested y-planes does'nt correspond parameters of projection for variable " + var + "(" + name + ")" + " (yb = " + ri->yb +
+                     ", ye = " + ri->ye + ", Ny = " + v.proj->Ny()};
+
+  if(ri->xe > nx - 1) return {pref, "Internal error: xe > nx - 1 for variable " + var + "(" + name + ")" + " (" + ri->xe + " > " + nx + " - 1)"};
+  if(ri->xb > nx - 1) return {pref, "Internal error: xb > nx - 1 for variable " + var + "(" + name + ")" + " (" + ri->xb + " > " + nx + " - 1)"};
+
+  const size_t rnx = (ri->xb <= ri->xe) ? (ri->xe - ri->xb + 1) : (nx + ri->xe - ri->xb + 1);
+  if(rnx != v.proj->Nx())
+   return {pref, "Internal error: number of requested x-planes does'nt correspond parameters of projection for variable " + var + "(" + name + ")" + " (xb = " + ri->xb +
+                     ", xe = " + ri->xe + ", Nx = " + v.proj->Nx()};
+ }
+
+ return {};
+}
+
+RetVal<std::shared_ptr<Data2D>> Adapter::Read2D(const MString& var, size_t it)
+{
+ static const MString pref = "Adapter::Read2D";
+
+ {
+  auto ret = ReadCheck(var, it);
+  if(!ret) return ret.Add(pref, "Check not passed");
+ }
+
+ // Clear cache
+ if(!cache2D) cache2D.reset(new decltype(cache2D)::element_type);
+ if(it != curtime) cache2D->clear();
+ curtime = it;
+
+ // Check cache
+ if(cache2D->contains(var)) return (*cache2D)[var];
+
+ const auto& v = vars[var];
+
+ if(!v.Read2D) return {pref, "No read function for the variable " + var};
+
+ auto ret = v.Read2D(*this, &v.info, v.proj, curit);
+ if(!ret) return ret.Add(pref, "Can't read variable " + var);
+ (*cache2D)[var] = ret.Value();
+ return (*cache2D)[var];
+}
+
+RetVal<std::shared_ptr<Data3D>> Adapter::Read3D(const MString& var, size_t it)
+{
+ static const MString pref = "Adapter::Read3D";
+
+ {
+  auto ret = ReadCheck(var, it);
+  if(!ret) return ret.Add(pref, "Check not passed");
+ }
+
+ // Clear cache
+ if(!cache3D) cache3D.reset(new decltype(cache3D)::element_type);
+ if(it != curtime) cache3D->clear();
+ curtime = it;
+
+ // Check cache
+ if(cache3D->contains(var)) return (*cache3D)[var];
+
+ const auto& v = vars[var];
+
+ if(!v.vert) return {pref, "Variable " + var + " is 2D in the this dataset"};
+ if(!v.Read3D) return {pref, "No read function for the variable " + var};
+
+ auto ret = v.Read3D(*this, &v.info, v.proj, v.vert, curit);
+ if(!ret) return ret.Add(pref, "Can't read variable " + var);
+ (*cache3D)[var] = ret.Value();
+ return (*cache3D)[var];
+}

src/CF.cpp

@@ -0,0 +1,469 @@
+#define MICHLIB_NOSOURCE
+#include "CF.h"
+
+using michlib::DetGeoDomain;
+
+MString CF::GetTVarName(const NCZarr& nc)
+{
+ for(const auto& v: nc.VarNames()) // Try to define coordinates by attribute standard_name or attribute axis
+ {
+  auto havestname = nc.HasAtt(v, "standard_name");
+  auto haveaxis   = nc.HasAtt(v, "axis");
+  if(!(havestname || haveaxis)) continue;
+
+  auto stname = nc.AttString(v, "standard_name");
+  auto axis   = nc.AttString(v, "axis");
+
+  if(stname == "time" || axis == "T") return v;
+ }
+
+ // If time not found just check variable "time"
+ if(nc.HasVar("time")) return "time";
+
+ return "";
+}
+
+struct CF::AxisInfo CF::GetAxisInfo(const NCZarr& nc)
+{
+ struct AxisInfo out;
+
+ for(const auto& v: nc.VarNames()) // Try to define coordinates by attribute standard_name or attribute axis
+ {
+  auto havestname = nc.HasAtt(v, "standard_name");
+  auto haveaxis   = nc.HasAtt(v, "axis");
+  if(!(havestname || haveaxis)) continue;
+
+  auto stname = nc.AttString(v, "standard_name");
+  auto axis   = nc.AttString(v, "axis");
+  bool islon = false, islat = false, isdepth = false, istime = false;
+  if(stname == "longitude") islon = true;
+  if(stname == "latitude") islat = true;
+  if(stname == "depth") isdepth = true;
+  if(stname == "time") istime = true;
+
+  if(!out.xvarname.Exist() && axis == "X") islon = true;
+  if(!out.yvarname.Exist() && axis == "Y") islat = true;
+  if(!out.zvarname.Exist() && axis == "Z") isdepth = true;
+  if(!out.tvarname.Exist() && axis == "T") istime = true;
+
+  if(islon) out.xvarname = v;
+  if(islat) out.yvarname = v;
+  if(isdepth) out.zvarname = v;
+  if(istime) out.tvarname = v;
+ }
+
+ // If time not found just check variable "time"
+ if(!out.tvarname.Exist() && nc.HasVar("time")) out.tvarname = "time";
+
+ // Dimensions
+ out.tdims = out.tvarname.Exist() ? nc.NDim(out.tvarname) : 0;
+ out.zdims = out.zvarname.Exist() ? nc.NDim(out.zvarname) : 0;
+ out.ydims = out.yvarname.Exist() ? nc.NDim(out.yvarname) : 0;
+ out.xdims = out.xvarname.Exist() ? nc.NDim(out.xvarname) : 0;
+
+ if(out.tdims == 1) out.tdimname = nc.DimNames(out.tvarname)[0];
+ if(out.zdims == 1) out.zdimname = nc.DimNames(out.zvarname)[0];
+ if(out.ydims == 1) out.ydimname = nc.DimNames(out.yvarname)[0];
+ if(out.xdims == 1) out.xdimname = nc.DimNames(out.xvarname)[0];
+
+ return out;
+}
+
+std::tuple<MDateTime, time_t, bool> CF::ParseRefDate(const MString& refdate)
+{
+ MDateTime out;
+ time_t    step = 0;
+
+ MString rstr;
+ auto    words = michlib::Split_on_words(refdate);
+ auto    ci    = words.begin();
+ if(ci != words.end())
+ {
+  if(*ci == "seconds") step = 1;
+  if(*ci == "minutes") step = 60;
+  if(*ci == "hours") step = 3600;
+  if(*ci == "days") step = 3600 * 24;
+  ci++;
+ }
+ if(ci != words.end()) ci++;              // skip "since"
+ if(ci != words.end()) rstr = *ci++;      // Day
+ if(ci != words.end()) rstr += " " + *ci; // Hours
+ bool success = out.FromString(rstr);
+ return {out, step, success};
+}
+
+RetVal<std::vector<MDateTime>> CF::ReadTimes(const NCZarr& nc, const MString& tname)
+{
+ static const MString   pref = "CF::ReadTimes";
+ std::vector<double>    time;
+ std::vector<MDateTime> out;
+
+ {
+  auto ret = nc.Read(tname, time);
+  if(!ret) return ret.Add(pref, "Can't read time");
+ }
+
+ MDateTime refdate;
+ time_t    step = 0;
+ {
+  auto units = nc.AttString(tname, "units");
+  if(!units.Exist()) return {pref, "Can't read refdate"};
+  auto [rd, st, suc] = ParseRefDate(units);
+  if(!suc) return {pref, "Can't parse " + units + " to refdate"};
+  if(st == 0) return {pref, "Can't get timestep from string " + units};
+  refdate = rd;
+  step    = st;
+ }
+
+ out.resize(time.size());
+ for(size_t i = 0; i < time.size(); i++) out[i] = refdate + static_cast<time_t>(time[i] * step);
+ return out;
+}
+
+Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad, std::unique_ptr<NCZarr>&& pnc)
+{
+ static const MString pref = "CF::FillAdapterFromCF";
+
+ bool debug = args.contains("debug");
+
+ auto axisinfo = GetAxisInfo(*pnc);
+
+ if(!axisinfo.tvarname.Exist()) return {pref, "Can't find time variable in the dataset"};
+ if(!axisinfo.yvarname.Exist()) return {pref, "Can't find latitude variable in the dataset"};
+ if(!axisinfo.xvarname.Exist()) return {pref, "Can't find longitude variable in the dataset"};
+
+ if(axisinfo.tdims != 1) return {pref, "Unsupported number of time dimensions: " + MString(axisinfo.tdims)};
+ if(axisinfo.ydims != 1) return {pref, "Unsupported number of latitude dimensions: " + MString(axisinfo.ydims)};
+ if(axisinfo.xdims != 1) return {pref, "Unsupported number of longitude dimensions: " + MString(axisinfo.xdims)};
+
+ if(axisinfo.zvarname.Exist() && axisinfo.zdims != 1) return {pref, "Unsupported number of depth dimensions: " + MString(axisinfo.zdims)};
+
+ if(debug)
+ {
+  pars.AddParameter("T variable", axisinfo.tvarname + "/" + axisinfo.tdims + "/" + axisinfo.tdimname);
+  pars.AddParameter("Z variable", axisinfo.zvarname + "/" + axisinfo.zdims + "/" + axisinfo.zdimname);
+  pars.AddParameter("Y variable", axisinfo.yvarname + "/" + axisinfo.ydims + "/" + axisinfo.ydimname);
+  pars.AddParameter("X variable", axisinfo.xvarname + "/" + axisinfo.xdims + "/" + axisinfo.xdimname);
+ }
+
+ std::shared_ptr<Projection> proj;
+ std::shared_ptr<Vertical>   vert;
+
+ // Forming ReadInfo
+ std::shared_ptr<Adapter::ReadInfo> rinfo(new Adapter::ReadInfo);
+ rinfo->xdname = axisinfo.xdimname;
+ rinfo->ydname = axisinfo.ydimname;
+ rinfo->zdname = axisinfo.zdimname;
+ rinfo->tdname = axisinfo.tdimname;
+
+ std::vector<double> depths;
+ if(axisinfo.zvarname.Exist())
+ {
+  auto ret = pnc->Read(axisinfo.zvarname, depths);
+  if(!ret) return ret.Add(pref, "Can't read depths");
+ }
+
+ // Filling xb, yb, xe, ye
+ // We assume that latitude and longitude are sorted in ascending order
+ {
+  std::vector<double> lons, lats;
+  {
+   auto ret = pnc->Read(axisinfo.xvarname, lons);
+   if(!ret) return ret.Add(pref, "Can't read longitudes");
+  }
+  {
+   auto ret = pnc->Read(axisinfo.yvarname, lats);
+   if(!ret) return ret.Add(pref, "Can't read latitudes");
+  }
+
+  real lonb, latb, lone, late;
+  bool global = lons.back() - lons.front() + 1.5 * (lons.back() - lons.front()) / (lons.size() - 1) > 360.0;
+  auto dom    = DetGeoDomain(lons.front(), lons.back());
+  lonb        = ToGeoDomain(args.contains("lonb") ? args.at("lonb").ToReal() : lons.front(), dom);
+  lone        = ToGeoDomain(args.contains("lone") ? args.at("lone").ToReal() : lons.back(), dom);
+  latb        = args.contains("latb") ? args.at("latb").ToReal() : lats.front();
+  late        = args.contains("late") ? args.at("late").ToReal() : lats.back();
+
+  if(global && lonb > lons.back()) return {pref, "Argument lonb must be lesser then " + MString(lons.back())};
+  if(global && lone < lons.front()) return {pref, "Argument lone must be greater then " + MString(lons.front())};
+  if(latb > lats.back()) return {pref, "Argument latb must be lesser then " + MString(lats.back())};
+  if(late < lats.front()) return {pref, "Argument late must be greater then " + MString(lats.front())};
+
+  if(global)
+  {
+   // Special case when the longitude lies in a small sector between the end and the start
+   if(lonb < lons.front()) lonb = lons.back();
+   if(lone > lons.back()) lone = lons.front();
+  }
+  else
+  {
+   if(lonb < lons.front()) lonb = lons.front();
+   if(lone > lons.back()) lone = lons.back();
+  }
+  if(latb < lats.front()) latb = lats.front();
+  if(late > lats.back()) late = lats.back();
+
+  for(rinfo->xb = 0; rinfo->xb + 1 < lons.size() && lonb >= lons[rinfo->xb + 1]; rinfo->xb++);
+  for(rinfo->yb = 0; rinfo->yb + 1 < lats.size() && latb >= lats[rinfo->yb + 1]; rinfo->yb++);
+
+  for(rinfo->xe = lons.size() - 1; rinfo->xe >= 1 && lone <= lons[rinfo->xe - 1]; rinfo->xe--);
+  for(rinfo->ye = lats.size() - 1; rinfo->ye >= 1 && late <= lats[rinfo->ye - 1]; rinfo->ye--);
+
+  if(debug)
+  {
+   pars.AddParameter("Requested lonb", lonb);
+   pars.AddParameter("Requested lone", lone);
+   pars.AddParameter("Requested latb", latb);
+   pars.AddParameter("Requested late", late);
+   pars.AddParameter("xb", rinfo->xb);
+   pars.AddParameter("xe", rinfo->xe);
+   pars.AddParameter("yb", rinfo->yb);
+   pars.AddParameter("ye", rinfo->ye);
+  }
+
+  if(rinfo->xb == rinfo->xe) return {pref, "Longitude interval too small"};
+  if(rinfo->yb == rinfo->ye) return {pref, "Latitude interval too small"};
+  if(rinfo->xe == 0) return {pref, "Sorry, special case then xe=0 does'nt handle. Increase lone."};
+  if(rinfo->xb == lons.size() - 1) return {pref, "Sorry, special case then xb=lons.size()-1 does'nt handle. Decrease lonb."};
+
+  lonb = lons[rinfo->xb];
+  lone = lons[rinfo->xe];
+  latb = lats[rinfo->yb];
+  late = lats[rinfo->ye];
+  if(lonb > lone)
+  {
+   if(lonb > 0.0)
+    lonb -= 360.0;
+   else
+    lone += 360.0;
+  }
+
+  pars.AddParameter("lonb", michlib::ToGeoDomainPos(lonb));
+  pars.AddParameter("lone", michlib::ToGeoDomainPos(lone));
+  pars.AddParameter("latb", latb);
+  pars.AddParameter("late", late);
+
+  decltype(lats) rlats = lats | std::views::drop(rinfo->yb) | std::views::take(rinfo->ye - rinfo->yb + 1) | std::ranges::to<std::vector>();
+  decltype(lons) rlons;
+
+  if(rinfo->xe >= rinfo->xb)
+   rlons = lons | std::views::drop(rinfo->xb) | std::views::take(rinfo->xe - rinfo->xb + 1) | std::ranges::to<std::vector>();
+  else
+  {
+   rlons.insert(rlons.begin(), lons.begin() + rinfo->xb, lons.end());
+   std::ranges::transform(lons.begin(), lons.begin() + rinfo->xe + 1, std::back_inserter(rlons), [](auto lon) { return lon + 360.0; });
+  }
+
+  // Check uniform and create projection
+  if(CF::CheckUniform(rlons, 0.0005) && CF::CheckUniform(rlats, 0.0005))
+  {
+   proj.reset(new Projection(Projection::Create<Projection::Type::EQC>(rlons, rlats)));
+   pars.AddParameter("lonstep", (rlons.back() - rlons.front()) / (rlons.size() - 1));
+   pars.AddParameter("latstep", (rlats.back() - rlats.front()) / (rlats.size() - 1));
+  }
+  else
+  {
+   michlib::errmessage("Warning! Nonuniform grid");
+   proj.reset(new Projection(Projection::Create<Projection::Type::IRREGULARXY>(std::move(rlons), std::move(rlats))));
+  }
+ }
+
+ // ReadInfo for the plain data
+ std::shared_ptr<Adapter::ReadInfo> rinfoplane(new Adapter::ReadInfo(*rinfo));
+ rinfoplane->zdname = "";
+ rinfoplane->zb = rinfoplane->ze = rinfo->zb = rinfo->ze = 0;
+
+ // Filling zb, ze
+ bool depthinv = false;
+ if(depths.size() > 1)
+ {
+  if(depths.back() <= 0 && depths.front() <= 0) std::ranges::transform(depths, depths.begin(), std::negate{});
+  if(depths.back() < depths.front() && depths.size() > 1)
+  {
+   depthinv = true;
+   for(size_t i = 0; i < depths.size() - i - 1; i++) std::swap(depths[i], depths[depths.size() - i - 1]);
+   rinfoplane->zb = rinfoplane->ze = rinfo->zb = rinfo->ze = depths.size() - 1;
+  }
+
+  if(args.contains("layer") || args.contains("depth"))
+  {
+   size_t layer = 0;
+   if(args.contains("layer"))
+    layer = args.at("layer").ToInteger<size_t>();
+   else
+   {
+    auto d = args.at("depth").ToReal();
+    if(d < depths.front()) d = depths.front();
+    if(d > depths.back()) d = depths.back();
+
+    for(size_t i = 0; i < depths.size() - 1; i++)
+     if(d >= depths[i] && d <= depths[i + 1])
+     {
+      layer = (d - depths[i] <= depths[i + 1] - d) ? i : (i + 1);
+      break;
+     }
+    if(debug) pars.AddParameter("Requested depth", d);
+   }
+   if(layer > depths.size() - 1) return {pref, "Layer must be lesser then " + MString(depths.size())};
+   rinfoplane->zb = rinfoplane->ze = rinfo->zb = rinfo->ze = depthinv ? depths.size() - layer - 1 : layer;
+
+   vert.reset(new Vertical(Vertical::Create<Vertical::Type::FIXED>(depths[layer])));
+
+   pars.AddParameter("layer", layer);
+   pars.AddParameter("depth", depths[layer]);
+   if(debug) pars.AddParameter("zb", rinfo->zb);
+  }
+  else if(args.contains("layers") || args.contains("depths"))
+  {
+   size_t layer1 = 0, layer2 = 0;
+
+   if(args.contains("layers"))
+   {
+    if(args.at("layers") == "ALL" || args.at("layers") == "all")
+    {
+     layer1 = 0;
+     layer2 = depths.size() - 1;
+    }
+    else
+    {
+     auto layerstr = args.at("layers").Split(":");
+
+     if(layerstr.size() > 2 || layerstr.size() == 0) return {pref, "Layers must have format layer1:layer2"};
+     if(layerstr.size() == 1)
+      layer1 = layer2 = layerstr[0].ToInteger<size_t>();
+     else
+     {
+      layer1 = layerstr[0].ToInteger<size_t>();
+      layer2 = layerstr[1].ToInteger<size_t>();
+     }
+    }
+   }
+   else
+   {
+    real d1 = 0.0, d2 = 0.0;
+    auto depthstr = args.at("depths").Split(":");
+    if(depthstr.size() > 2 || depthstr.size() == 0) return {pref, "Depths must have format depth1:depth2"};
+    if(depthstr.size() == 1)
+     d1 = d2 = depthstr[0].ToReal();
+    else
+    {
+     d1 = depthstr[0].ToReal();
+     d2 = depthstr[1].ToReal();
+    }
+    if(d1 < depths.front()) d1 = depths.front();
+    if(d1 > depths.back()) d1 = depths.back();
+    if(d2 < depths.front()) d2 = depths.front();
+    if(d2 > depths.back()) d2 = depths.back();
+
+    for(size_t i = 0; i < depths.size() - 1; i++)
+     if(d1 >= depths[i] && d1 <= depths[i + 1])
+     {
+      layer1 = (d1 - depths[i] <= depths[i + 1] - d1) ? i : (i + 1);
+      break;
+     }
+    for(size_t i = 0; i < depths.size() - 1; i++)
+     if(d2 >= depths[i] && d2 <= depths[i + 1])
+     {
+      layer2 = (d2 - depths[i] <= depths[i + 1] - d2) ? i : (i + 1);
+      break;
+     }
+    if(debug)
+    {
+     pars.AddParameter("Requested depth1", d1);
+     pars.AddParameter("Requested depth2", d2);
+    }
+   }
+
+   if(layer1 > layer2) return {pref, "Layer2 must be greater then layer1"};
+   if(layer1 > depths.size() - 1) return {pref, "Layer1 must be lesser then " + MString(depths.size())};
+   if(layer2 > depths.size() - 1) return {pref, "Layer2 must be lesser then " + MString(depths.size())};
+   rinfoplane->zb = rinfoplane->ze = depthinv ? depths.size() - layer1 - 1 : layer1;
+   rinfo->zb                       = depthinv ? depths.size() - layer1 - 1 : layer1;
+   rinfo->ze                       = depthinv ? depths.size() - layer2 - 1 : layer2;
+
+   pars.AddParameter("layer1", layer1);
+   pars.AddParameter("depth1", depths[layer1]);
+   pars.AddParameter("layer2", layer2);
+   pars.AddParameter("depth2", depths[layer2]);
+
+   if(debug)
+   {
+    pars.AddParameter("zb", rinfo->zb);
+    pars.AddParameter("ze", rinfo->ze);
+   }
+
+   // Create vertical projection
+   {
+    decltype(depths) rdepths(depths.begin() + layer1, depths.begin() + layer2 + 1);
+    vert.reset(new Vertical(Vertical::Create<Vertical::Type::HORIZONTS>(rdepths)));
+   }
+  }
+  else
+  {
+   pars.AddParameter("layer", 0);
+   pars.AddParameter("depth", depths[0]);
+   if(debug) pars.AddParameter("zb", rinfo->zb);
+
+   vert.reset(new Vertical(Vertical::Create<Vertical::Type::FIXED>(depths[0])));
+  }
+ }
+ if(debug) pars.AddParameter("depthinv", depthinv);
+
+ //michlib::message("depth1=", depths[depthinv ? depths.size() - 1 - rinfo->zb : rinfo->zb], ", zb=", rinfo->zb);
+ //michlib::message("depth2=", depths[depthinv ? depths.size() - 1 - rinfo->ze : rinfo->ze], ", ze=", rinfo->ze);
+ //michlib::message(depthinv ? "Depths inverted" : "Depths not inverted");
+
+ // Read variables information
+ for(const auto& v: pnc->VarNames())
+ {
+  Adapter::VInfo vinfo;
+  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;
+
+  vinfo.name = v;
+  vinfo.unit = "1";
+  if(dnames.size() == 4)
+   vinfo.rinfo = rinfo;
+  else
+   vinfo.rinfo = rinfoplane;
+
+  if(pnc->HasAtt(v, "units")) vinfo.unit = pnc->AttString(v, "units");
+  if(pnc->HasAtt(v, "standard_name")) vinfo.stname = pnc->AttString(v, "standard_name");
+  if(pnc->HasAtt(v, "long_name")) vinfo.lname = pnc->AttString(v, "long_name");
+  vinfo.offset  = pnc->HasAtt(v, "add_offset") ? pnc->AttReal(v, "add_offset") : 0.0;
+  vinfo.scale   = pnc->HasAtt(v, "scale_factor") ? pnc->AttReal(v, "scale_factor") : 1.0;
+  vinfo.fillval = std::visit(
+      [](auto r) -> real
+      {
+       if constexpr(std::is_convertible_v<decltype(r), real>)
+        return static_cast<real>(r);
+       else
+        return NAN;
+      },
+      pnc->VarFill(v));
+
+  vinfo.targetunit = vinfo.unit;
+  if(vinfo.stname.SubStr(vinfo.stname.Len() - 7, 8) == "velocity") vinfo.targetunit = "cm/s";
+
+  MString name = StName2Name(vinfo.stname).Exist() ? StName2Name(vinfo.stname) : v;
+  if(dnames.size() == 3) { ad.Add2DVariable(name, std::move(vinfo), proj); }
+  else { ad.Add3DVariable(name, std::move(vinfo), proj, vert); }
+ }
+
+ return Error();
+}
+
+/*
+RetVal<std::shared_ptr<Data2D>> CF::PTemp3DReader(const Adapter& ad, const void* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it)
+{
+ static const MString pref = "CF::PTemp3DReader";
+
+ auto v = michlib::pointer_cast<const struct Adapter::VInfo*>(vinfo);
+
+ const struct Adapter::ReadInfo* ri = v->rinfo.get();
+
+ std::shared_ptr<Data3D> out(new Data3D(proj, vert, {.unit = v->targetunit, .stname = v->stname, .lname = v->lname, .comment = v->comment}));
+}
+*/

src/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(EXENAME odm)
 set(ACTIONLISTINC ${CMAKE_CURRENT_BINARY_DIR}/../include/actionlist.h) # Include actions files and define the actions classes list
 set(DATALISTINC ${CMAKE_CURRENT_BINARY_DIR}/../include/datalist.h) # Include data sources files and define the data sources classes list
+set(USAGEH ${CMAKE_CURRENT_BINARY_DIR}/../include/usage.h)
 
 find_package(PkgConfig REQUIRED)
 find_library(netcdf netcdf REQUIRED)
@@ -21,10 +22,18 @@ file(GLOB srcs CONFIGURE_DEPENDS *.cpp)
 add_executable(${EXENAME} ${srcs} ${ACTIONLISTINC} ${SOURCELISTINC})
 target_include_directories(${EXENAME} PRIVATE ../michlib/michlib ${CMAKE_CURRENT_BINARY_DIR}/../include)
 
-target_link_libraries(${EXENAME} ${linker_options} ${netcdf} ${udunits} OpenMP::OpenMP_CXX CURL::libcurl ${JSONCPP_LINK_LIBRARIES} ${BLOSC_LINK_LIBRARIES} ${LIBPQ_LINK_LIBRARIES} LibXml2::LibXml2 SQLite::SQLite3 teos)
+target_link_libraries(${EXENAME} PRIVATE ${linker_options} ${netcdf} ${udunits} OpenMP::OpenMP_CXX CURL::libcurl ${JSONCPP_LINK_LIBRARIES} ${BLOSC_LINK_LIBRARIES} ${LIBPQ_LINK_LIBRARIES} LibXml2::LibXml2 SQLite::SQLite3 teos)
 set_target_properties(${EXENAME} PROPERTIES POSITION_INDEPENDENT_CODE ON)
 install(TARGETS ${EXENAME})
 
+add_custom_command(OUTPUT ${USAGEH}
+ COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/../genusageh ${CMAKE_CURRENT_SOURCE_DIR}/../doc.txt ${USAGEH}
+ DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/../doc.txt ${CMAKE_CURRENT_SOURCE_DIR}/../genusageh
+ VERBATIM
+)
+add_custom_target(usageh DEPENDS ${USAGEH})
+add_dependencies(${EXENAME} usageh)
+
 # Begin generation of actions list
 get_target_property(INCS ${EXENAME} INCLUDE_DIRECTORIES)
 list(TRANSFORM INCS PREPEND "-I")

src/conscolors.cpp

@@ -0,0 +1,4 @@
+#define MICHLIB_NOSOURCE
+#include "conscolors.h"
+
+bool C::tty = false;

src/copcat.cpp

@@ -77,6 +77,25 @@ RetVal<MString> CopernicusCatalog::ProductURL(const MString& prod) const
  return {pref, "unknown product: " + prod};
 }
 
+RetVal<MString> CopernicusCatalog::ProductTittle(const MString& prod) const
+{
+ static const MString pref = "CopernicusCatalog::ProductTtitle";
+
+ MString url;
+ {
+  auto ret = ProductURL(prod);
+  if(!ret) return ret.Add(pref, "Can't get url for the product " + prod);
+  url = ret.Value();
+ }
+
+ auto ret = GetJSON(url);
+ if(!ret) return ret.Add(pref, "Can't download product " + prod);
+
+ const auto& title = ret.Value()["title"];
+ if(title.type() != Json::stringValue) return {pref, "no valid \"title\" value in the product " + prod + " description"};
+ return MString(title.asString().c_str());
+}
+
 RetVal<std::vector<MString>> CopernicusCatalog::DatasetList(const MString& prod) const
 {
  static const MString pref = "CopernicusCatalog::DatasetList";
@@ -135,6 +154,25 @@ RetVal<MString> CopernicusCatalog::DatasetURL(const MString& prod, const MString
  return {pref, "unknown dataset: " + dataset};
 }
 
+RetVal<MString> CopernicusCatalog::DatasetTitle(const MString& prod, const MString& dataset) const
+{
+ static const MString pref = "CopernicusCatalog::DatasetTitle";
+
+ MString url;
+ {
+  auto ret = DatasetURL(prod, dataset);
+  if(!ret) return ret.Add(pref, "Can't get url for the dataset " + dataset);
+  url = ret.Value();
+ }
+
+ auto ret = GetJSON(url);
+ if(!ret) return ret.Add(pref, "Can't download dataset " + dataset);
+
+ const auto& title = ret.Value()["properties"]["title"];
+ if(!title || title.type() != Json::stringValue) return {pref, "title for the dataset " + dataset + " not found"};
+ return MString(title.asString().c_str());
+}
+
 RetVal<MString> CopernicusCatalog::AssetURL(const MString& prod, const MString& dataset, const MString& asset) const
 {
  static const MString pref = "CopernicusCatalog::AssetURL";

src/ncfilew.cpp

@@ -141,3 +141,227 @@ MString NCFileW::AddVariable(const MString& name, const MString& stname, const M
 
  return "";
 }
+
+michlib::Error NCFileW::AddTimeData(const std::vector<MDateTime>& times)
+{
+ static const MString pref = "NCFileW::AddTimeData";
+
+ static constexpr auto                                            stepunits = std::to_array<uint>({3600 * 24, 3600, 60, 1});
+ static constexpr std::array<const std::string, stepunits.size()> stepnames{"days", "hours", "minutes", "seconds"};
+
+ MDateTime                  refdate;
+ size_t                     unitind(stepunits.size());
+ std::vector<michlib::int4> steps(times.size());
+
+ if(steps.size() == 0) return {pref, "No time data available"};
+ refdate = times[0];
+ unitind = 0;
+
+ for(size_t i = 0; i < steps.size(); i++)
+ {
+  auto delta = times[i] - refdate;
+  while(delta.Seconds() % stepunits[unitind] != 0) unitind++;
+ }
+
+ for(size_t i = 0; i < steps.size(); i++)
+ {
+  auto delta = times[i] - refdate;
+  steps[i]   = michlib::int_cast<decltype(steps)::value_type>(delta.Seconds() / stepunits[unitind]);
+ }
+
+ MString unit = stepnames[unitind].c_str() + (" since " + refdate.ToString());
+
+ AddDim("time", steps.size());
+ AddVar("time", Type2NCType<decltype(steps)::value_type>, "time");
+ SetComp("time", compress);
+ AddAtt("time", "standard_name", "time");
+ AddAtt("time", "long_name", "time");
+ AddAtt("time", "units", unit);
+ WriteVar("time", steps.data());
+
+ if(!*this) return {pref, MString("Can't add time data to the netcdf file: ") + ErrMessage()};
+
+ return {};
+}
+
+michlib::Error NCFileW::Create(const MString& name, const Adapter& ad, const michlib_internal::ParameterListEx& pars, int compression)
+{
+ static const MString pref = "NCFileW::Create";
+ if(Ok()) return {pref, "File already opened"};
+
+ const float node_offset = 0.0;
+
+ Open(name);
+ if(!*this) return {pref, "Can't create netcdf file " + name + ": " + ErrMessage()};
+
+ MString rets;
+
+ rets = AddAtt("node_offset", node_offset);
+ if(rets.Exist()) return {pref, rets};
+
+ rets = AddAtts(pars);
+ if(rets.Exist()) return {pref, rets};
+
+ compress = compression;
+
+ {
+  auto ret = AddTimeData(ad.IndexedTimes());
+  if(!ret) return ret.Add(pref, "Can't write time data");
+ }
+
+ for(const auto& v: ad.Vars())
+ {
+  bool is2D = v.second.Read2D ? true : false;
+  bool is3D = v.second.Read3D ? true : false;
+
+  const auto& vname   = v.first;                     // Variable name
+  const auto& proj    = *v.second.proj;              // Horizontal projection
+  const auto& vert    = *v.second.vert;              // Vertical projection
+  const auto& units   = v.second.info.targetunit;    // Unit
+  const auto& stname  = v.second.info.stname;        // Standard name
+  const auto& lname   = v.second.info.lname;         // Long name
+  const auto& comment = v.second.info.comment;       // Comment
+  const auto& xdname  = v.second.info.rinfo->xdname; // X-dimension name
+  const auto& ydname  = v.second.info.rinfo->ydname; // Y-dimension name
+  const auto& zdname  = v.second.info.rinfo->zdname; // Z-dimension name
+  const auto& fillval = v.second.info.fillval;       // _FillValue
+
+  if((is2D && is3D) || (!is2D && !is3D)) return {pref, "Can't determine type of variable " + vname};
+
+  // Z-dimension
+  if(is3D && vert.IsValid())
+  {
+   // TODO Поддержка неизотропных сеток
+   if(!vert.IsIsotropic()) return {pref, "Non-isotropic vertical grids still unsupported"};
+   // Check if z-dimension already exist
+   int  did;
+   auto st = nc_inq_dimid(*this, zdname.Buf(), &did);
+   if(st != NC_NOERR) // Add zdata
+   {
+    AddDim(zdname, vert.Nz());
+    AddVar(zdname, Type2NCType<float>, zdname);
+    SetComp(zdname, compress);
+    AddAtt(zdname, "standard_name", "depth");
+    AddAtt(zdname, "long_name", "Depth");
+    AddAtt(zdname, "units", "m");
+    std::vector<float> data(vert.Nz());
+    for(size_t i = 0; i < vert.Nz(); i++) data[i] = static_cast<float>(vert.Depth(i));
+    WriteVar(zdname, data.data());
+    if(!*this) return {pref, MString("Can't add depth data to the netcdf file: ") + ErrMessage()};
+   }
+   else // Check compatibility with zdata in file
+   {
+    size_t len;
+    st = nc_inq_dimlen(*this, did, &len);
+    if(st != NC_NOERR) return {pref, "Can't get size of depth dimension"};
+    if(len != vert.Nz()) return {pref, "Length of depth dimension does'nt corrrespond to number of vertical layers of variable " + vname};
+   }
+  }
+
+  // X- and Y-dimensions
+  // TODO Поддержка непрямоугольных сеток
+  if(!proj.IsRectangular()) return {pref, "Non-rectangular horizontal grids still unsupported"};
+  {
+   int  lonid, latid;
+   auto stlon = nc_inq_dimid(*this, xdname.Buf(), &lonid);
+   auto stlat = nc_inq_dimid(*this, ydname.Buf(), &latid);
+
+   if(stlon != NC_NOERR) // Add longitudes
+   {
+    AddDim(xdname, proj.Nx());
+    AddVar(xdname, Type2NCType<float>, xdname);
+    SetComp(xdname, compress);
+    AddAtt(xdname, "standard_name", "longitude");
+    AddAtt(xdname, "long_name", "Longitude");
+    AddAtt(xdname, "units", "degrees_east");
+    std::vector<float> data(proj.Nx());
+    for(size_t i = 0; i < proj.Nx(); i++) data[i] = static_cast<float>(proj.Lon(i));
+    WriteVar(xdname, data.data());
+    if(!*this) return {pref, MString("Can't add longitude data to the netcdf file: ") + ErrMessage()};
+   }
+   else // Check compatibility with longitudes in file
+   {
+    size_t len;
+    auto   st = nc_inq_dimlen(*this, lonid, &len);
+    if(st != NC_NOERR) return {pref, "Can't get size of longitude dimension"};
+    if(len != proj.Nx()) return {pref, "Length of longitude dimension does'nt corrrespond to number of longitudes for variable " + vname};
+   }
+
+   if(stlat != NC_NOERR) // Add latitudes
+   {
+    AddDim(ydname, proj.Ny());
+    AddVar(ydname, Type2NCType<float>, ydname);
+    SetComp(ydname, compress);
+    AddAtt(ydname, "standard_name", "latitude");
+    AddAtt(ydname, "long_name", "Latitude");
+    AddAtt(ydname, "units", "degrees_north");
+    std::vector<float> data(proj.Ny());
+    for(size_t i = 0; i < proj.Ny(); i++) data[i] = static_cast<float>(proj.Lat(i));
+    WriteVar(ydname, data.data());
+    if(!*this) return {pref, MString("Can't add latitude data to the netcdf file: ") + ErrMessage()};
+   }
+   else // Check compatibility with longitudes in file
+   {
+    size_t len;
+    auto   st = nc_inq_dimlen(*this, latid, &len);
+    if(st != NC_NOERR) return {pref, "Can't get size of latitude dimension"};
+    if(len != proj.Ny()) return {pref, "Length of latitude dimension does'nt corrrespond to number of latitudes for variable " + vname};
+   }
+  }
+
+  // Variable
+  if(is2D) AddVar(vname, Type2NCType<float>, "time", ydname, xdname);
+  if(is3D) AddVar(vname, Type2NCType<float>, "time", zdname, ydname, xdname);
+  SetComp(vname, compress);
+  if(stname.Exist()) AddAtt(vname, "standard_name", stname);
+  if(lname.Exist()) AddAtt(vname, "long_name", lname);
+  if(units.Exist()) AddAtt(vname, "units", units);
+  if(comment.Exist()) AddAtt(vname, "comment", comment);
+  AddAtt(vname, "_FillValue", static_cast<float>(fillval));
+  if(!*this) return {pref, MString("Can't add variable " + vname + " to the netcdf file: ") + ErrMessage()};
+ } // end vars for
+
+ return {};
+}
+
+template<class D>
+requires(std::is_same_v<D, Data2D> || std::is_same_v<D, Data3D>)
+michlib::Error NCFileW::WriteVariable(const MString& name, const D& data, size_t it)
+{
+ static constexpr const bool   is2D = std::is_same_v<D, Data2D>;
+ static const MString          pref = "NCFileW::WriteVariable";
+ static constexpr const size_t ND   = is2D ? 3 : 4;
+
+ // Check variable presence
+ int varid;
+ {
+  auto ret = nc_inq_varid(*this, name.Buf(), &varid);
+  if(ret != NC_NOERR) return {pref, "Can't find variable " + name};
+ }
+
+ // Check is variable 2D
+ {
+  int  ndim;
+  auto ret = nc_inq_var(*this, varid, nullptr, nullptr, &ndim, nullptr, nullptr);
+  if(ret != NC_NOERR) return {pref, "Can't get number of dimensions for the variable " + name};
+  if(ndim != ND) return {pref, "Number of dimensions for the variable " + name + " is " + ndim + ", but must be " + ND};
+ }
+
+ // Check dimensions?
+
+ // Write
+ const double* pdata;
+ if constexpr(is2D)
+  pdata = &data(0, 0);
+ else
+  pdata = &data(0, 0, 0);
+
+ //WriteVar(name, it, michlib::pointer_cast<const float*>(pdata));
+ WriteVar(name, it, pdata);
+ if(!*this) return {pref, MString("Failed to write variable " + name + ": ") + ErrMessage()};
+
+ return {};
+}
+
+template michlib::Error NCFileW::WriteVariable<Data2D>(const MString& name, const Data2D& data, size_t it);
+template michlib::Error NCFileW::WriteVariable<Data3D>(const MString& name, const Data3D& data, size_t it);

src/odm.cpp

@@ -1,57 +1,13 @@
 #include "odm.h"
 
+const char usage[] =
+#include "usage.h"
+    ;
+
 inline void Usage(const MString& arg0)
 {
  message(arg0 + " (<key>=<value>...)");
- message("Keys are:");
- message(" action. What the program should do. May be: info, tsc, uv, genintfile. Default: info.");
- message(" Keys for action=info. Print some information about dataset.");
- message("  source. Required. May be: NEMO, NEMOBIO, HYCOM, AVISO, AVISOLOCAL, BINFILE");
- message("  Keys for source=NEMO");
- message("   dataset. Can be DT, NRT or NRT6. Default: DT");
- message("  Keys for source=NEMOBIO");
- message("   dataset. Can be DT or NRT. Default: DT");
- message("  Keys for source=HYCOM");
- message("   dataset. Can be Forecast, Hindcast or Reanalysis. Default: Forecast");
- message("  Keys for source=AVISO");
- message("   dataset. Can be DT, NRT, EckmanDT or EckmanNRT. Default: DT");
- message(" Keys for action=tsc. Get temperature, salinity, chlorofill from dataset.");
- message("  source. Required. May be: NEMO, NEMOBIO, HYCOM, AVISO, AVISOLOCAL, BINFILE");
- message("  var. Required. May be: U, U2, u, v, temp, ptemp, pdens, sal, chl, mld, ssh, w, NO3, PO4, Si, prprod, Cchl.");
- message("  time. Time moment or regular expression. If present, timeb and timee must be absent");
- message("  timeb, timee. Time interval. If present, time must be absent");
- message("   lonb, lone, latb, late. Required. Region of interest");
- message("  out. Output file. Default: out.bin");
- message("  Keys for source=NEMO");
- message("   dataset. Can be DT, NRT or NRT6. Default: DT");
- message("   layer and/or depth. Layer or depth of NEMO dataset. If depth is specified, layer is ignored. Both ignored, if var=mld or var=ssh. Default: layer=0");
- message("  Keys for source=NEMOBIO");
- message("   dataset. Can be DT or NRT. Default: DT");
- message("   layer and/or depth. Layer or depth of NEMOBIO dataset. If depth is specified, layer is ignored. Default: layer=0");
- message("  Keys for source=HYCOM");
- message("   dataset. Can be Forecast, Hindcast or Reanalysis. Default: Forecast");
- message("   layer and/or depth. Layer or depth of HYCOM dataset. If depth is specified, layer is ignored. Both ignored, if var=mld or var=ssh. Default: layer=0");
- message("  Keys for source=AVISO");
- message("   dataset. Can be DT, NRT, EckmanDT or EckmanNRT. Default: DT");
- message("   layer and/or depth. Layer or depth of AVISO dataset. If depth is specified, layer is ignored. Both ignored for datasets DT and NRT. Default: layer=0");
- message("  Keys for source=BINFILE");
- message("   dataset. Path or DataID of interpolation file");
- message(" Keys for action=uv. Get velocity field and its derivatives.");
- message("  source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL, BINFILE");
- message("  time. Time moment or regular expression. If present, timeb and timee must be absent");
- message("  timeb, timee. Time interval. If present, time must be absent");
- message("  out. Output file for components of velocity field, divergency, rotor and Okubo-Weiss parameter. If absent, this data not calculated.");
- message("  velout. Output file for sparse velocity field. If present, parameters shiftx, shifty, skipx, skipy controls the sparsing.");
- message("  shiftx, shifty. Shift of the sparsed grid. Used only if velout parameter is present. Default: 0");
- message("  skipx, skipy. How many poinst skipped in the sparsed grid. Used only if velout parameter is present. Default: 1, output each point");
- message("  stpout. Output file for stationary points. If absent, this data not calculated.");
- message("  Keys for specific sources see in the section action=tsc.");
- message(" Keys for action=genintfile. Create file with velocity field for interpolation.");
- message("  source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL");
- message("  time. Time moment or regular expression. If present, timeb and timee must be absent");
- message("  timeb, timee. Time interval. If present, time must be absent");
- message("  out. Output file, required.");
- message("  Keys for specific sources see in the section action=tsc.");
+ message(usage);
 }
 
 int main(int argc, char** argv)
@@ -69,6 +25,11 @@ int main(int argc, char** argv)
 
  auto args = ParseArgs(argc, argv);
 
+ if((args.contains("color") && !args.contains("nocolor")) || (!args.contains("color") && args.contains("nocolor")))
+  C::Init(args.contains("color"));
+ else
+  C::Init();
+
  MString action = args.contains("action") ? args["action"] : "info";
 
  Data data;

src/udunitscpp.cpp

@@ -0,0 +1,5 @@
+#define MICHLIB_NOSOURCE
+#include "udunitscpp.h"
+
+ut_system* UtUnit::system = nullptr;
+size_t     UtUnit::count  = 0;