Compare commits

...

15 Commits

  1. 9
      CMakeLists.txt
  2. 4
      actions/actiongenintfile.h
  3. 6
      actions/actionint.h
  4. 7
      actions/actionmirror.h
  5. 4
      include/actiondep.h
  6. 135
      include/cache.h
  7. 54
      include/copcat.h
  8. 9
      include/curlfuncs.h
  9. 2
      include/layereddata.h
  10. 207
      include/layereddataz.h
  11. 10
      include/mirrorfuncs.h
  12. 10
      include/ncfuncs.h
  13. 87
      include/ncsimple.h
  14. 836
      include/nczarrcommon.h
  15. 175
      include/zarr.h
  16. 2
      michlib
  17. 6
      sources/AVISO.h
  18. 2
      sources/AVISOLOCAL.h
  19. 2
      sources/BINFILE.h
  20. 146
      sources/COPERNICUS.cpp
  21. 21
      sources/COPERNICUS.h
  22. 20
      sources/NEMO.h
  23. 6
      sources/NEMOBIO.h
  24. 10
      sources/VYLET.cpp
  25. 6
      src/CMakeLists.txt
  26. 168
      src/copcat.cpp
  27. 298
      src/layereddataz.cpp
  28. 60
      src/mirrorfuncs.cpp
  29. 104
      src/ncfuncs.cpp
  30. 231
      src/ncsimple.cpp
  31. 249
      src/zarr.cpp

9
CMakeLists.txt

@ -39,7 +39,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_INTERPROCEDURAL_OPTIMIZATION ON)
include_directories(include sources sources-add GSW-C)
add_compile_options(-Wall)
add_compile_options(-Wall -Wno-deprecated-declarations)
# Dwarf-4 support check
include(CheckCXXCompilerFlag)
@ -54,6 +54,13 @@ if(COMPILER_SUPPORTS_ANALYZER AND STATIC_ANALYZER)
add_compile_options(-fanalyzer)
endif()
# Disable michlib warnings
add_compile_options(-Wno-deprecated-declarations)
CHECK_CXX_COMPILER_FLAG(-Wno-class-memaccess COMPILER_SUPPORTS_CLASSMEMACCESS)
if(COMPILER_SUPPORTS_CLASSMEMACCESS)
add_compile_options(-Wno-class-memaccess)
endif()
add_library(teos STATIC GSW-C/gsw_oceanographic_toolbox.c GSW-C/gsw_saar.c)
set_target_properties(teos PROPERTIES LINKER_LANGUAGE C)

4
actions/actiongenintfile.h

@ -69,7 +69,7 @@ template<class D> MString ActionGenIntFile::DoAction(const CLArgs& args, D& ds)
fw.SetParameter("dx", data.XStep() * 60.0);
fw.SetParameter("dy", data.YStep() * 60.0);
fw.SetParameter("nt", tindexes.size());
fw.SetParameter("dt", (ds.Time(tindexes[1]) - ds.Time(tindexes[0])) / 86400.0);
fw.SetParameter("dt", (ds.Time(tindexes[1]) - ds.Time(tindexes[0])).D());
fw.SetParameter("FillValue", data.Fillval());
fw.UsePrefix("Info");
fw.SetParameter("Mode", mode ? "from ssh" : "normal");
@ -92,7 +92,7 @@ template<class D> MString ActionGenIntFile::DoAction(const CLArgs& args, D& ds)
fw.SetParameter("dotxdotylonlat2v", "(1.852/0.864)*%doty");
fw.SetParameter("BeginDate", ds.Time(tindexes.front()).ToTString());
fw.SetParameter("EndDate", ds.Time(tindexes.back()).ToTString());
fw.SetParameter("Timestep", ds.Time(tindexes[1]) - ds.Time(tindexes[0]));
fw.SetParameter("Timestep", (ds.Time(tindexes[1]) - ds.Time(tindexes[0])).Seconds());
fw.SetParameter("DataID", michlib::UniqueId());
fw.SetParameter("Creation command", args.at("_cmdline"));
}

6
actions/actionint.h

@ -187,9 +187,9 @@ template<class D> MString ActionINT::DoAction(const CLArgs& args, D& ds)
for(size_t j = 0; j < temp.size(); j++) hi.emplace_back(std::move(temp[j]));
}
}
auto step = ds.Time(hiind) - ds.Time(loind);
auto delta = points[i].t - ds.Time(loind);
real trel = static_cast<real>(delta) / static_cast<real>(step);
auto step = (ds.Time(hiind) - ds.Time(loind)).S();
auto delta = (points[i].t - ds.Time(loind)).S();
real trel = delta / step;
fw.Write(points[i].lon0);
fw.Write(points[i].lat0);

7
actions/actionmirror.h

@ -3,21 +3,20 @@
#include "merrors.h"
using michlib::message;
using michlib::Error;
template<class T>
concept MirrorSupported = requires(T t, const CLArgs& args) {
{
t.Mirror(args)
} -> std::convertible_to<MString>;
} -> std::convertible_to<Error>;
};
ADD_ACTION(Mirror, mirror, MirrorSupported<Source>);
template<class D> MString ActionMirror::DoAction(const CLArgs& args, D& data)
{
//auto resop = data.Open(args);
//if(resop.Exist()) return "Can't open source: " + resop;
auto res = data.Mirror(args);
if(res.Exist()) return "Mirroring failed: " + res;
if(!res) return "Mirroring failed";
return "";
};

4
include/actiondep.h

@ -41,13 +41,13 @@ struct TimeData
for(size_t i = 0; i < steps.size(); i++)
{
auto delta = data.Time(tindexes[i]) - refdate;
while(delta % stepunits[unitind] != 0) unitind++;
while(delta.Seconds() % stepunits[unitind] != 0) unitind++;
}
for(size_t i = 0; i < steps.size(); i++)
{
auto delta = data.Time(tindexes[i]) - refdate;
steps[i] = michlib::int_cast<decltype(steps)::value_type>(delta / stepunits[unitind]);
steps[i] = michlib::int_cast<decltype(steps)::value_type>(delta.Seconds() / stepunits[unitind]);
}
}

135
include/cache.h

@ -1,10 +1,13 @@
#pragma once
#include "MString.h"
#include "merrors.h"
#include <libpq-fe.h>
#include <sqlite3.h>
#include <time.h>
#include <variant>
using michlib::int_cast;
using michlib::MString;
using michlib::pointer_cast;
class GenericCache
{
@ -121,6 +124,129 @@ class SQLiteCache: public GenericCache
explicit operator bool() const { return db != nullptr; }
};
class PostgreSQLCache: public GenericCache
{
PGconn* conn = nullptr;
bool CheckCon() const
{
if(!*this) return false;
if(PQstatus(conn) == CONNECTION_OK) return true;
PQreset(conn);
return PQstatus(conn) == CONNECTION_OK;
}
template<class D> static D Invert(D d)
{
using michlib::int1;
if(sizeof(D) <= 1) return d;
D out;
int1* pout = pointer_cast<int1*>(&out);
int1* pin = pointer_cast<int1*>(&d);
for(size_t i = 0; i < sizeof(D); i++) pout[sizeof(D) - i - 1] = pin[i];
return out;
}
public:
bool Init(const MString& name)
{
Close();
conn = PQconnectdb(name.Buf());
if(PQstatus(conn) != CONNECTION_OK)
{
michlib::errmessage(PQerrorMessage(conn));
Close();
return false;
}
// Create table
if(false)
{
auto* res = PQexec(conn, "CREATE TABLE IF NOT EXISTS cache(key TEXT PRIMARY KEY NOT NULL, value BYTEA, exptime TIMESTAMP(0) NOT NULL);");
if(PQresultStatus(res) != PGRES_COMMAND_OK)
{
michlib::errmessage(PQresStatus(PQresultStatus(res)));
michlib::errmessage(PQerrorMessage(conn));
PQclear(res);
Close();
}
else
PQclear(res);
}
return true;
}
void Close()
{
if(conn != nullptr) PQfinish(conn);
conn = nullptr;
}
virtual bool Put(const MString& key, const MString& value, size_t ttl) const override
{
if(!CheckCon()) return false;
auto interval = michlib::int_cast<michlib::int8>(ttl);
michlib::int8 rinterval = Invert(interval);
const char* params[] = {key.Buf(), value.Buf(), pointer_cast<const char*>(&rinterval)};
int plens[] = {int_cast<int>(key.Len()), int_cast<int>(value.Len()), 8};
int pfor[] = {0, 1, 1};
auto* res = PQexecParams(conn,
"INSERT INTO cache(key,value,exptime) VALUES($1,$2,localtimestamp + ($3::bigint ||' seconds')::interval)"
"ON CONFLICT(key) DO UPDATE SET value=EXCLUDED.value, exptime=EXCLUDED.exptime;",
3, nullptr, params, plens, pfor, 1);
if(PQresultStatus(res) != PGRES_COMMAND_OK)
{
michlib::errmessage(PQresStatus(PQresultStatus(res)));
michlib::errmessage(PQerrorMessage(conn));
PQclear(res);
return false;
}
PQclear(res);
return true;
}
virtual std::pair<MString, bool> Get(const MString& key) const override
{
if(!CheckCon()) return {"", false};
const char* params[] = {key.Buf()};
int plens[] = {int_cast<int>(key.Len())};
int pfor[] = {0};
auto* res = PQexecParams(conn, "SELECT value from cache WHERE key=$1::text AND exptime>localtimestamp;", 1, nullptr, params, plens, pfor, 1);
if(PQresultStatus(res) != PGRES_TUPLES_OK)
{
michlib::errmessage(PQresStatus(PQresultStatus(res)));
michlib::errmessage(PQerrorMessage(conn));
PQclear(res);
return {"", false};
}
else if(PQntuples(res) == 0)
{
PQclear(res);
return {"", false};
}
MString val(PQgetvalue(res, 0, 0), PQgetlength(res, 0, 0));
PQclear(res);
return {std::move(val), true};
}
virtual ~PostgreSQLCache() override
{
if(!CheckCon()) return;
auto* res = PQexec(conn, "DELETE FROM cache WHERE exptime<localtimestamp;");
PQclear(res);
Close();
}
explicit operator bool() const { return conn != nullptr; }
};
inline GenericCache* CreateCache(const MString& cachedesc)
{
auto i = cachedesc.GetPos(':');
@ -140,6 +266,13 @@ inline GenericCache* CreateCache(const MString& cachedesc)
if(*ret) return ret;
delete ret;
}
if(name == "postgre" || name == "postgres" || name == "postgresql")
{
auto ret = new PostgreSQLCache;
ret->Init(par);
if(*ret) return ret;
delete ret;
}
return nullptr;
}

54
include/copcat.h

@ -0,0 +1,54 @@
#pragma once
#include "cache.h"
#include "curlfuncs.h"
#include "merrors.h"
#include <json/json.h>
using michlib::Error;
using michlib::RetVal;
class CopernicusCatalog
{
static const MString caturl;
std::unique_ptr<GenericCache> cache;
CURLRAII chandle;
Json::Value catalog;
// Download catalog
Error GetCatalog();
// Asset url from dataset
RetVal<MString> AssetURL(const MString& prod, const MString& dataset, const MString& asset) const;
public:
CopernicusCatalog();
// Download JSON from url
RetVal<Json::Value> GetJSON(const MString& url) const;
// List of products
RetVal<std::vector<MString>> ProductList() const;
// List of datasets in product
RetVal<std::vector<MString>> DatasetList(const MString& prod) const;
// URL of product
RetVal<MString> ProductURL(const MString& prod) const;
// URL of dataset
RetVal<MString> DatasetURL(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"); }
// URL of timechuncked data (files) in dataset
RetVal<MString> DatasetTimeURL(const MString& prod, const MString& dataset) const { return AssetURL(prod, dataset, "timeChunked"); }
// URL of geochuncked data (files) in dataset
RetVal<MString> DatasetGeoURL(const MString& prod, const MString& dataset) const { return AssetURL(prod, dataset, "geoChunked"); }
bool Valid() const { return catalog.isObject(); }
explicit operator bool() const { return Valid(); }
};

9
include/curlfuncs.h

@ -14,9 +14,16 @@ class CURLRAIIDT
class CURLRAII: public std::unique_ptr<CURL, CURLRAIIDT>
{
char err[CURL_ERROR_SIZE];
public:
CURLRAII() { reset(curl_easy_init()); }
CURLRAII()
{
reset(curl_easy_init());
curl_easy_setopt(*this, CURLOPT_ERRORBUFFER, err);
}
operator CURL*() const { return get(); }
const char* Err() const { return err; }
};
// Curl writeback function, write to MString

2
include/layereddata.h

@ -162,7 +162,7 @@ class LayeredData: public NCFuncs
return times[i];
}
time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
time_t Timestep() const { return isOk() ? (times[1] - times[0]).Seconds() : 0; }
MString Title() const { return title; }

207
include/layereddataz.h

@ -0,0 +1,207 @@
#pragma once
#include "gsw.h"
#include "ncfuncs.h"
#include "simple2ddata.h"
#include "zarr.h"
#include <memory>
using michlib::Ceil;
using michlib::DetGeoDomain;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;
class LayeredDataZ: public NCFuncs
{
public:
using Data = Simple2DData;
private:
class NC: public Zarr
{
std::vector<MDateTime> times;
public:
Error ReadTimes(const MString& tname)
{
static const MString pref = "LayeredDataZ::NC::ReadTimes";
if(!*this) return Error(pref, "Dataset not open");
std::vector<double> time;
{
auto ret = Read(tname, time);
if(!ret) return ret.Add(pref, "Can't read time");
}
MDateTime refdate;
time_t step = 0;
{
auto units = AttString(tname, "units");
if(!units.Exist()) return Error(pref, "Can't read refdate");
auto [rd, st, suc] = Refdate(units);
if(!suc) return Error(pref, "Can't parse " + units + " to refdate");
if(st == 0) return Error(pref, "Can't get timestep from string " + units);
refdate = rd;
step = st;
}
times.resize(time.size());
for(size_t i = 0; i < time.size(); i++) times[i] = refdate + static_cast<time_t>(time[i] * step);
return Error();
}
MDateTime Begin() const { return times.front(); }
MDateTime End() const { return times.back(); }
const std::vector<MDateTime>& Times() const { return times; }
size_t Index(MDateTime tm) const
{
if(tm < Begin() || tm > End()) return 0;
size_t b = 0, e = times.size() - 1;
if(tm == times[b]) return b + 1;
if(tm == times[e]) return e + 1;
while(e - b > 1)
{
size_t c = (e + b) / 2;
if(tm == times[c]) return c + 1;
if(tm > times[c])
b = c;
else
e = c;
}
return 0;
}
};
std::vector<NC> nc;
std::vector<real> depths;
bool depthinv;
std::vector<MDateTime> times;
struct CoordNames dname;
real lonb, latb, lone, late;
real lonstep, latstep;
MString title;
class EnvVar
{
MString name, oldvalue;
bool activated, saved;
public:
EnvVar(): activated(false) {}
~EnvVar() { Deactivate(); }
void Activate(const MString& var, const MString& val)
{
if(activated) Deactivate();
name = var;
char* curval = getenv(name.Buf());
if(nullptr == curval)
saved = false;
else
{
oldvalue = curval;
saved = true;
}
setenv(name.Buf(), val.Buf(), 1);
}
void Deactivate()
{
if(!activated) return;
if(saved)
setenv(name.Buf(), oldvalue.Buf(), 1);
else
unsetenv(name.Buf());
activated = false;
}
};
EnvVar proxy;
protected:
struct Parameters: public BaseParameters
{
size_t xb, yb, xe, ye, layer;
virtual ~Parameters() override = default;
};
// TODO: RetVal
MString Open(const MString& dataset);
void SetTitle(const MString& newtitle) { title = newtitle; }
public:
MString Info() const;
std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const;
bool Read(const MString& vname, std::map<MString, Data>& cache, const BaseParameters* ip, size_t i) const;
bool isOk() const { return nc.size() > 0; }
explicit operator bool() const { return nc.size() > 0; }
real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; }
real Depth(const BaseParameters* ip) const { return Depth(dynamic_cast<const struct Parameters*>(ip)->layer); }
real Lon(size_t ix) const { return isOk() ? (lonb + ix * lonstep) : -1000.0; }
real Lat(size_t iy) const { return isOk() ? (latb + iy * latstep) : -1000.0; }
size_t NDepths() const { return depths.size(); }
size_t NTimes() const { return times.size(); }
MDateTime Time(size_t i) const
{
if(!isOk() || i >= times.size()) return MDateTime();
return times[i];
}
time_t Timestep() const { return isOk() ? (times[1] - times[0]).Seconds() : 0; }
MString Title() const { return title; }
MString Dump(const struct Parameters* ppar) const
{
// clang-format off
return
"Current settings:\n" + MString() +
" Longitudes: from " + Lon(ppar->xb) + " (" + ppar->xb + ") to "+ Lon(ppar->xe) + " (" + ppar->xe + ")\n" +
" Latitudes: from " + Lat(ppar->yb) + " (" + ppar->yb + ") to "+ Lat(ppar->ye) + " (" + ppar->ye + ")\n" +
" Depth: layer " + ppar->layer + ", depth " + Depth(ppar->layer) + " m\n";
// clang-format on
}
VarPresence CheckVar(const MString& vname) const
{
return NCFuncs::CheckVar(vname, [this](const MString& vn) { return HaveVar(vn); });
}
private:
template<class DataType> Data ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct Parameters* p) const;
bool HaveVar(const MString& vname) const
{
for(size_t i = 0; i < nc.size(); i++)
if(NCFuncs::HaveVar(nc[i], vname)) return true;
return false;
}
std::tuple<MString, size_t, size_t> VarNameLoc(const MString vname, MDateTime tm) const
{
for(size_t i = 0; i < nc.size(); i++)
{
auto tind = nc[i].Index(tm);
if(tind == 0) continue;
for(const auto& v: nc[i].Vars())
{
auto stname = nc[i].AttString(v.Name(), "standard_name");
if(!stname.Exist()) continue;
if(StName2Name(stname) == vname) return {v.Name(), i, tind - 1};
}
}
return {"", 0, 0};
}
};

10
include/mirrorfuncs.h

@ -8,6 +8,8 @@
#include <vector>
using michlib::MDateTime;
using michlib::RetVal;
using michlib::Error;
class DIRRAIIDT
{
@ -50,13 +52,13 @@ inline MString FileName(const MString& name)
bool MakePath(const MString& dname);
// Get local file list
std::pair<std::vector<struct FileInfo>, MString> ReadLocalFileList(const MString& dir, const MString& path = "");
RetVal<std::vector<struct FileInfo>> ReadLocalFileList(const MString& dir, const MString& path = "");
// Download file to the local mirror
MString DownloadFile(const CURLRAII& chandle, const struct FileInfo& rinfo, const MString& root);
Error DownloadFile(const CURLRAII& chandle, const struct FileInfo& rinfo, const MString& root);
// Remove file from the local mirror
MString RemoveFile(const struct FileInfo& linfo);
Error RemoveFile(const struct FileInfo& linfo);
// Updare file in the local mirror
MString UpdateFile(const CURLRAII& chandle, const struct FileInfo& rinfo, const struct FileInfo& linfo, const MString& root);
Error UpdateFile(const CURLRAII& chandle, const struct FileInfo& rinfo, const struct FileInfo& linfo, const MString& root);

10
include/ncfuncs.h

@ -2,6 +2,8 @@
#include "DataAdapters/ncfilealt.h"
#include "basedata.h"
#include "mdatetime.h"
#include "ncsimple.h"
#include "zarr.h"
#include <map>
#include <set>
#include <tuple>
@ -27,7 +29,13 @@ class NCFuncs
static CoordNames GetCNames(const NCFileA& nc);
static CoordNames GetDNames(const NCFileA& nc);
static bool HaveVar(const NCFileA& nc, const MString& vname);
template<class HV> static VarPresence CheckVar(const MString& vname, HV hv)
template<class NcZarrFunctions> static void GetVars(const NcZarrRead<NcZarrFunctions>& nc, std::set<MString>& vars);
template<class NcZarrFunctions> static CoordNames GetCNames(const NcZarrRead<NcZarrFunctions>& nc);
template<class NcZarrFunctions> static CoordNames GetDNames(const NcZarrRead<NcZarrFunctions>& nc);
template<class NcZarrFunctions> static bool HaveVar(const NcZarrRead<NcZarrFunctions>& nc, const MString& vname);
template<class HV> static VarPresence CheckVar(const MString& vname, HV hv)
{
if(!hv(vname))
{

87
include/ncsimple.h

@ -0,0 +1,87 @@
#pragma once
#include "nczarrcommon.h"
#include <netcdf.h>
#include <variant>
class NCSimpleTypes: public NcZarrTypes
{
protected:
template<class VType> class ReadedData
{
using Vec = std::vector<size_t>;
private:
std::unique_ptr<VType[]> data;
public:
ReadedData() = default;
ReadedData(std::unique_ptr<VType[]>&& d): data(std::move(d)) {}
VType operator()(size_t lini) const { return data[lini]; }
};
};
class NCSimpleFunctions: public NCSimpleTypes
{
int ncid;
RetVal<std::vector<Attribute>> ReadAtts(int vid) const;
protected:
NCSimpleFunctions(): ncid(0) {}
template<class VType> RetVal<ReadedData<VType>> Read(const MString& var, const size_t* start, const size_t* count) const
{
static const MString pref = "NCSimpleFunctions::Read";
size_t ind = FindInd(var, vars);
const size_t N = vars[ind].NDim();
std::unique_ptr<VType[]> cdata;
size_t dsize = 1;
for(size_t i = 0; i < N; i++) dsize *= count[i];
cdata.reset(new VType[dsize]);
int vid;
int res = nc_inq_varid(ncid, var.Buf(), &vid);
if(res != NC_NOERR) return Error(pref, MString("nc_inq_varid error: ") + nc_strerror(res));
if constexpr(std::is_same_v<VType, float>)
res = nc_get_vara_float(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, double>)
res = nc_get_vara_double(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, int>)
res = nc_get_vara_int(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, long>)
res = nc_get_vara_long(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, long long>)
res = nc_get_vara_longlong(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, short>)
res = nc_get_vara_short(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, signed char>)
res = nc_get_vara_schar(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, unsigned int>)
res = nc_get_vara_uint(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, unsigned long long>)
res = nc_get_vara_ulonglong(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, unsigned short>)
res = nc_get_vara_ushort(ncid, vid, start, count, cdata.get());
else if constexpr(std::is_same_v<VType, unsigned char>)
res = nc_get_vara_ubyte(ncid, vid, start, count, cdata.get());
else
return Error(pref, "Unsupported variable type");
if(res != NC_NOERR) return Error(pref, MString("nc_get_vara error: ") + nc_strerror(res));
return ReadedData<VType>(std::move(cdata));
}
public:
~NCSimpleFunctions() { nc_close(ncid); }
Error Open(const MString& filename);
};
using NCSimple = NcZarrRead<NCSimpleFunctions>;

836
include/nczarrcommon.h

@ -0,0 +1,836 @@
#pragma once
#include "merrors.h"
#include <utility>
#include <variant>
using michlib::Error;
using michlib::int1;
using michlib::int2;
using michlib::int4;
using michlib::int8;
using michlib::int_cast;
using michlib::MString;
using michlib::RetVal;
using michlib::uint1;
using michlib::uint8;
class NcZarrTypes
{
protected:
using AttVT = std::variant<std::monostate, int8, uint8, double, MString, bool>;
class ArrCounter
{
using VT = std::vector<size_t>;
const VT count;
VT ind;
bool end;
public:
static size_t Index(const VT& i, const VT& c)
{
size_t out = 0;
size_t mul = 1;
for(size_t ii = i.size(); ii != 0; ii--)
{
out += mul * i[ii - 1];
mul *= c[ii - 1];
}
return out;
}
static VT Index(size_t lind, const VT& c)
{
VT out(c.size());
size_t j = lind;
for(auto i = c.size(); i > 0; i--)
{
out[i - 1] = j % c[i - 1];
j = j / c[i - 1];
}
return out;
}
ArrCounter() = delete;
ArrCounter(const VT& cnt): count(cnt), ind(cnt.size(), 0), end(false) {}
size_t operator[](size_t i) const { return ind[i]; }
ArrCounter& operator++()
{
size_t curind = count.size();
while(curind != 0)
{
ind[curind - 1]++;
if(ind[curind - 1] >= count[curind - 1])
{
ind[curind - 1] = 0;
curind--;
}
else
return *this;
}
ind = count;
end = true;
return *this;
}
explicit operator bool() const { return !end; }
size_t Index() const { return Index(ind, count); }
size_t Index(const VT& i) const { return Index(i, count); }
VT Index(size_t lind) const { return Index(lind, count); }
size_t Count(size_t i) const { return count[i]; }
const VT& VIndex() const { return ind; }
VT VIndex(const VT& start) const
{
VT out(ind.size());
for(size_t i = 0; i < ind.size(); i++) out[i] = ind[i] + start[i];
return out;
}
const auto& Count() const { return count; }
size_t N() const
{
size_t out = 1;
for(size_t i = 0; i < count.size(); i++) out *= count[i];
return out;
}
};
public:
enum class AttType
{
UNDEF,
INT,
UINT,
REAL,
STRING,
BOOL
};
enum class VarType
{
UNDEF,
FLOAT,
DOUBLE,
INT1,
INT2,
INT4,
INT8,
UINT1
};
protected:
template<VarType VT, class Dummy = void> struct VarType2Type;
template<class Dummy> struct VarType2Type<VarType::FLOAT, Dummy>
{
using type = float;
};
template<class Dummy> struct VarType2Type<VarType::DOUBLE, Dummy>
{
using type = double;
};
template<class Dummy> struct VarType2Type<VarType::INT1, Dummy>
{
using type = int1;
};
template<class Dummy> struct VarType2Type<VarType::INT2, Dummy>
{
using type = int2;
};
template<class Dummy> struct VarType2Type<VarType::INT4, Dummy>
{
using type = int4;
};
template<class Dummy> struct VarType2Type<VarType::INT8, Dummy>
{
using type = int8;
};
template<class Dummy> struct VarType2Type<VarType::UINT1, Dummy>
{
using type = uint1;
};
template<VarType VT> using Type = VarType2Type<VT>::type;
static constexpr size_t SizeOf(VarType vt)
{
switch(vt)
{
case(VarType::UNDEF): return 0;
case(VarType::FLOAT): return sizeof(Type<VarType::FLOAT>);
case(VarType::DOUBLE): return sizeof(Type<VarType::DOUBLE>);
case(VarType::INT1): return sizeof(Type<VarType::INT1>);
case(VarType::INT2): return sizeof(Type<VarType::INT2>);
case(VarType::INT4): return sizeof(Type<VarType::INT4>);
case(VarType::INT8): return sizeof(Type<VarType::INT8>);
case(VarType::UINT1): return sizeof(Type<VarType::UINT1>);
}
return 0;
}
template<class T> static size_t FindInd(const MString& name, const std::vector<T>& arr)
{
for(size_t i = 0; i < arr.size(); i++)
if(arr[i].Name() == name) return i;
return arr.size();
}
class Attribute: public AttVT
{
MString name;
public:
Attribute(const MString& n, AttVT&& v): AttVT(std::move(v)), name(n) {}
Attribute(const std::string& n, AttVT&& v): AttVT(std::move(v)), name(n.c_str(), n.size()) {}
const MString& Name() const { return name; }
AttType Type() const
{
if(std::holds_alternative<int8>(*this))
return AttType::INT;
else if(std::holds_alternative<uint8>(*this))
return AttType::UINT;
else if(std::holds_alternative<double>(*this))
return AttType::REAL;
else if(std::holds_alternative<MString>(*this))
return AttType::STRING;
else if(std::holds_alternative<bool>(*this))
return AttType::BOOL;
return AttType::UNDEF;
}
int8 I() const
{
if(std::holds_alternative<int8>(*this))
return std::get<int8>(*this);
else if(std::holds_alternative<uint8>(*this))
return int_cast<int8>(std::get<uint8>(*this));
else if(std::holds_alternative<double>(*this))
return static_cast<int8>(std::get<double>(*this));
else if(std::holds_alternative<MString>(*this))
return std::get<MString>(*this).ToInteger<int8>();
else if(std::holds_alternative<bool>(*this))
return std::get<bool>(*this) ? 1 : 0;
return 0;
}
uint8 U() const
{
if(std::holds_alternative<int8>(*this))
return int_cast<uint8>(std::get<int8>(*this));
else if(std::holds_alternative<uint8>(*this))
return std::get<uint8>(*this);
else if(std::holds_alternative<double>(*this))
return static_cast<uint8>(std::get<double>(*this));
else if(std::holds_alternative<MString>(*this))
return std::get<MString>(*this).ToInteger<uint8>();
else if(std::holds_alternative<bool>(*this))
return std::get<bool>(*this) ? 1 : 0;
return 0;
}
double D() const
{
if(std::holds_alternative<int8>(*this))
return std::get<int8>(*this);
else if(std::holds_alternative<uint8>(*this))
return std::get<uint8>(*this);
else if(std::holds_alternative<double>(*this))
return std::get<double>(*this);
else if(std::holds_alternative<MString>(*this))
return michlib_internal::RealType<sizeof(double)>::String2Real(std::get<MString>(*this).Buf());
else if(std::holds_alternative<bool>(*this))
return std::get<bool>(*this) ? 1 : 0;
return 0;
}
MString S() const
{
if(std::holds_alternative<int8>(*this))
return MString().FromInt(std::get<int8>(*this));
else if(std::holds_alternative<uint8>(*this))
return MString().FromUInt(std::get<uint8>(*this));
else if(std::holds_alternative<double>(*this))
return MString().FromReal(std::get<double>(*this));
else if(std::holds_alternative<MString>(*this))
return std::get<MString>(*this);
else if(std::holds_alternative<bool>(*this))
return MString().FromBool(std::get<bool>(*this));
return "";
}
bool B() const
{
if(std::holds_alternative<int8>(*this))
return std::get<int8>(*this) != 0;
else if(std::holds_alternative<uint8>(*this))
return std::get<uint8>(*this) != 0;
else if(std::holds_alternative<double>(*this))
return std::get<double>(*this) != 0.0;
else if(std::holds_alternative<MString>(*this))
return std::get<MString>(*this).ToBool();
else if(std::holds_alternative<bool>(*this))
return std::get<bool>(*this);
return false;
}
};
class Dimension
{
MString name;
size_t size;
public:
Dimension(const MString& str, size_t num): name(str), size(num) {}
const MString& Name() const { return name; }
size_t Size() const { return size; }
};
class Variable
{
public:
using FillType = std::variant<std::monostate, int8, uint8, double>;
private:
MString name;
VarType type = VarType::UNDEF;
std::vector<size_t> dims;
std::vector<Attribute> atts;
FillType fill;
public:
Variable(const MString& name_, VarType type_, std::vector<size_t>&& dims_, std::vector<Attribute>&& atts_, FillType fill_ = 0):
name(name_), type(type_), dims(std::move(dims_)), atts(std::move(atts_)), fill(fill_)
{
}
explicit operator bool() const { return type != VarType::UNDEF; }
const auto& Dims() const { return dims; }
size_t NDim() const { return dims.size(); }
size_t NAtt() const { return atts.size(); }
auto AttNames() const
{
std::vector<MString> out;
std::transform(atts.cbegin(), atts.cend(), std::back_inserter(out), [](const Attribute& a) { return a.Name(); });
return out;
}
AttType AttT(const MString& name) const
{
size_t ind = FindInd(name, atts);
return ind < atts.size() ? atts[ind].Type() : AttType::UNDEF;
}
int8 AttInt(const MString& name) const
{
size_t ind = FindInd(name, atts);
return ind < atts.size() ? atts[ind].I() : 0;
}
uint8 AttUInt(const MString& name) const
{
size_t ind = FindInd(name, atts);
return ind < atts.size() ? atts[ind].U() : 0;
}
double AttReal(const MString& name) const
{
size_t ind = FindInd(name, atts);
return ind < atts.size() ? atts[ind].D() : 0.0;
}
MString AttString(const MString& name) const
{
size_t ind = FindInd(name, atts);
return ind < atts.size() ? atts[ind].S() : MString();
}
bool AttBool(const MString& name) const
{
size_t ind = FindInd(name, atts);
return ind < atts.size() ? atts[ind].B() : false;
}
const MString& Name() const { return name; }
auto Type() const { return type; }
const auto& Fill() const { return fill; }
};
protected:
std::vector<Attribute> gats;
std::vector<Dimension> dims;
std::vector<Variable> vars;
public:
operator bool() const { return !vars.empty(); }
size_t NDim() const { return dims.size(); }
size_t NDim(const MString& var) const
{
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].NDim() : 0;
}
size_t NAtt() const { return gats.size(); }
auto AttNames() const
{
std::vector<MString> out;
std::transform(gats.cbegin(), gats.cend(), std::back_inserter(out), [](const Attribute& a) { return a.Name(); });
return out;
}
size_t NAtt(const MString& var) const
{
if(!var.Exist()) return NAtt();
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].NAtt() : 0;
}
auto AttNames(const MString& var) const
{
if(!var.Exist()) return AttNames();
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].AttNames() : decltype(AttNames())();
}
auto VarNames() const
{
std::vector<MString> out;
std::transform(vars.cbegin(), vars.cend(), std::back_inserter(out), [](const Variable& v) { return v.Name(); });
return out;
}
VarType VarT(const MString& var) const
{
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].Type() : VarType::UNDEF;
}
auto VarFill(const MString& var) const
{
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].Fill() : Variable::FillType();
}
auto DimNames() const
{
std::vector<MString> out;
std::transform(dims.cbegin(), dims.cend(), std::back_inserter(out), [](const Dimension& d) { return d.Name(); });
return out;
}
auto DimNames(const MString& var) const
{
size_t ind = FindInd(var, vars);
std::vector<MString> out;
if(ind >= vars.size()) return out;
auto vdims = vars[ind].Dims();
std::transform(vdims.cbegin(), vdims.cend(), std::back_inserter(out), [&dims = std::as_const(dims)](const size_t& i) { return dims[i].Name(); });
return out;
}
size_t DimSize(const MString& dim) const
{
size_t ind = FindInd(dim, dims);
return ind < dims.size() ? dims[ind].Size() : 0;
}
AttType AttT(const MString& var, const MString& name) const
{
if(!var.Exist())
{
size_t ind = FindInd(name, gats);
return ind < gats.size() ? gats[ind].Type() : AttType::UNDEF;
}
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].AttT(name) : AttType::UNDEF;
}
int8 AttInt(const MString& var, const MString& name) const
{
if(!var.Exist())
{
size_t ind = FindInd(name, gats);
return ind < gats.size() ? gats[ind].I() : 0;
}
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].AttInt(name) : 0;
}
uint8 AttUInt(const MString& var, const MString& name) const
{
if(!var.Exist())
{
size_t ind = FindInd(name, gats);
return ind < gats.size() ? gats[ind].U() : 0;
}
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].AttUInt(name) : 0;
}
double AttReal(const MString& var, const MString& name) const
{
if(!var.Exist())
{
size_t ind = FindInd(name, gats);
return ind < gats.size() ? gats[ind].D() : 0.0;
}
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].AttReal(name) : 0.0;
}
MString AttString(const MString& var, const MString& name) const
{
if(!var.Exist())
{
size_t ind = FindInd(name, gats);
return ind < gats.size() ? gats[ind].S() : MString();
}
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].AttString(name) : MString();
}
bool AttBool(const MString& var, const MString& name) const
{
if(!var.Exist())
{
size_t ind = FindInd(name, gats);
return ind < gats.size() ? gats[ind].B() : false;
}
size_t ind = FindInd(var, vars);
return ind < vars.size() ? vars[ind].AttBool(name) : false;
}
auto AttT(const MString& name) const { return AttT("", name); }
auto AttInt(const MString& name) const { return AttInt("", name); }
auto AttUInt(const MString& name) const { return AttUInt("", name); }
auto AttReal(const MString& name) const { return AttReal("", name); }
auto AttString(const MString& name) const { return AttString("", name); }
auto AttBool(const MString& name) const { return AttBool("", name); }
bool HasDim(const MString& name) const { return FindInd(name, dims) < dims.size(); }
bool HasVar(const MString& name) const { return FindInd(name, vars) < vars.size(); }
bool HasAtt(const MString& vname, const MString& aname) const { return AttT(vname, aname) != AttType::UNDEF; }
bool HasAtt(const MString& aname) const { return AttT(aname) != AttType::UNDEF; }
const auto& Vars() const { return vars; }
const auto& Dims() const { return dims; }
};
class DimReqDef
{
protected:
struct DimReq
{
static const auto fill = std::numeric_limits<size_t>::max();
MString name;
size_t beg, count;
DimReq(): name(MString()), beg(fill), count(fill) {}
DimReq(const char* n): name(n), beg(fill), count(fill) {}
DimReq(const MString& n): name(n), beg(fill), count(fill) {}
DimReq(MString&& n): name(std::move(n)), beg(fill), count(fill) {}
DimReq(const char* n, size_t s): name(n), beg(s), count(fill) {}
DimReq(const MString& n, size_t s): name(n), beg(s), count(fill) {}
DimReq(MString&& n, size_t s): name(std::move(n)), beg(s), count(fill) {}
DimReq(const char* n, size_t s, size_t c): name(n), beg(s), count(c) {}
DimReq(const MString& n, size_t s, size_t c): name(n), beg(s), count(c) {}
DimReq(MString&& n, size_t s, size_t c): name(std::move(n)), beg(s), count(c) {}
const MString& Name() const { return name; }
};
};
template<class C> class NcZarrRead: public C, public DimReqDef
{
template<class Data> static constexpr size_t Dimensionity()
{
if constexpr(requires(Data& d) { d(0, 0, 0, 0); }) return 4;
if constexpr(requires(Data& d) { d(0, 0, 0); }) return 3;
if constexpr(requires(Data& d) { d(0, 0); }) return 2;
if constexpr(requires(Data& d) { d(0); }) return 1;
return 0;
}
template<class Data, size_t D, class Dummy = void> struct DataTypeExtractorS;
template<class Data, class Dummy> struct DataTypeExtractorS<Data, 1, Dummy>
{
using type = std::decay_t<decltype(std::declval<Data>()(0))>;
};
template<class Data, class Dummy> struct DataTypeExtractorS<Data, 2, Dummy>
{
using type = std::decay_t<decltype(std::declval<Data>()(0, 0))>;
};
template<class Data, class Dummy> struct DataTypeExtractorS<Data, 3, Dummy>
{
using type = std::decay_t<decltype(std::declval<Data>()(0, 0, 0))>;
};
template<class Data, class Dummy> struct DataTypeExtractorS<Data, 4, Dummy>
{
using type = std::decay_t<decltype(std::declval<Data>()(0, 0, 0, 0))>;
};
template<class Data> using DataTypeExtractor = DataTypeExtractorS<Data, Dimensionity<Data>()>::type;
template<class VType, class Data, class Transform>
Error Read(const MString& vname, const std::vector<size_t>& transindex, Data& data, Transform transform, std::vector<DimReq> reqs) const
{
size_t nval = 1;
for(const auto& r: reqs) nval *= r.count;
const size_t indim = reqs.size();
constexpr size_t outdim = Dimensionity<Data>();
std::vector<size_t> start;
std::vector<size_t> count;
start.resize(indim);
count.resize(indim);
for(size_t i = 0; i < indim; i++)
{
start[i] = reqs[i].beg;
count[i] = reqs[i].count;
}
using DataType = DataTypeExtractor<Data>;
DataType fillout;
bool havefill = C::VarFill(vname).index() > 0;
VType fillin = std::visit(
[](auto v)
{
if constexpr(std::is_convertible_v<decltype(v), VType>)
return static_cast<VType>(v);
else
return std::numeric_limits<VType>::max();
},
C::VarFill(vname));
if constexpr(requires(Data& d) { // Data have own fillvalue
{
d.Fillval()
} -> std::convertible_to<DataType>;
})
fillout = data.Fillval();
else // Data does'nt have own fillvalue, using variable fillvalue
fillout = static_cast<DataType>(fillin);
auto ret = C::template Read<VType>(vname, start.data(), count.data());
if(!ret) return ret;
const auto& rawdata = ret.Value();
std::vector<size_t> mul(indim, 1);
for(size_t i = indim - 1; i > 0; i--) mul[i - 1] = mul[i] * count[i];
size_t inind = 0;
for(typename C::ArrCounter i(count); i; ++i)
{
// TODO: Remove this testing block
size_t cind = 0;
for(size_t j = 0; j < indim; j++) cind += i[j] * mul[j];
if(cind != inind) return {"NcZarrRead::Read", "Internal error"};
if(i.Index() != inind) return {"NcZarrRead::Read", "Internal error"};
if(inind != i.Index(i.Index(inind, count), count)) return {"NcZarrRead::Read", "Internal error"};
DataType out;
const VType& in = rawdata(inind);
if(havefill && in == fillin)
out = fillout;
else
out = transform(in);
if constexpr(outdim == 1)
data(i[transindex[0]]) = out;
else if constexpr(outdim == 2)
data(i[transindex[0]], i[transindex[1]]) = out;
else if constexpr(outdim == 3)
data(i[transindex[0]], i[transindex[1]], i[transindex[2]]) = out;
else if constexpr(outdim == 4)
data(i[transindex[0]], i[transindex[1]], i[transindex[2]], i[transindex[3]]) = out;
inind++;
}
return Error();
}
public:
// Request is string
template<class Data, class Transform> Error Read(const MString& vname, Data& data, Transform transform, const char* request) const
{
return Read(vname, data, transform, MString(request));
}
// Request by one dimension
template<class Data, class Transform> Error Read(const MString& vname, Data& data, Transform transform, DimReq&& req1) const
{
return Read(vname, data, transform, std::vector<DimReq>{std::move(req1)});
}
// Request by two dimension
template<class Data, class Transform> Error Read(const MString& vname, Data& data, Transform transform, DimReq&& req1, DimReq&& req2) const
{
return Read(vname, data, transform, std::vector<DimReq>{std::move(req1), std::move(req2)});
}
// Request by three dimension
template<class Data, class Transform> Error Read(const MString& vname, Data& data, Transform transform, DimReq&& req1, DimReq&& req2, DimReq&& req3) const
{
return Read(vname, data, transform, std::vector<DimReq>{std::move(req1), std::move(req2), std::move(req3)});
}
// Request by four dimension
template<class Data, class Transform> Error Read(const MString& vname, Data& data, Transform transform, DimReq&& req1, DimReq&& req2, DimReq&& req3, DimReq&& req4) const
{
return Read(vname, data, transform, std::vector<DimReq>{std::move(req1), std::move(req2), std::move(req3), std::move(req4)});
}
// Request full variable
template<class Data, class Transform> Error Read(const MString& vname, Data& data, Transform transform) const
{
static const MString pref = "NcZarrRead::Read";
if(!C::HasVar(vname)) return {pref, "Variable " + vname + " not found"};
std::vector<struct DimReq> pdims;
const auto vdims = C::DimNames(vname);
std::transform(
vdims.cbegin(), vdims.cend(), std::back_inserter(pdims), [this](const MString& n) -> struct DimReq {
return {n, 0, C::DimSize(n)};
});
return Read(vname, data, transform, pdims);
}
// Base function for all Read's
template<class Data, class Transform> Error Read(const MString& vname, Data& data, Transform transform, std::vector<DimReq> reqs) const
{
static const MString pref = "NcZarrRead::Read";
if(!C::HasVar(vname)) return {pref, "Variable " + vname + " not found"};
std::vector<struct DimReq> pdims;
{
const auto vdims = C::DimNames(vname);
std::transform(
vdims.cbegin(), vdims.cend(), std::back_inserter(pdims), [](const MString& n) -> struct DimReq {
return {n, 0, 1};
});
}
std::vector<size_t> transindex;
// Parse request
if(reqs.size() == 0) return {pref, "Empty request"};
for(const auto& req: reqs)
{
size_t ind = C::FindInd(req.name, pdims);
if(ind >= pdims.size()) return {pref, "Variable " + vname + " has no dimension " + req.name};
for(size_t i = 0; i < transindex.size(); i++)
if(transindex[i] == ind) return {pref, "Parameters for dimension " + req.name + " already defined"};
transindex.push_back(ind);
size_t dlen = C::DimSize(pdims[ind].name);
if(req.beg == req.fill && req.count == req.fill) // Only name, so, we request full length
{
pdims[ind].beg = 0;
pdims[ind].count = dlen;
}
else if(req.count == req.fill) // Name and first index
{
pdims[ind].beg = req.beg;
pdims[ind].count = 1;
}
else // Name, first index, count
{
pdims[ind].beg = req.beg;
pdims[ind].count = req.count;
}
// Sanity checks
if(pdims[ind].count <= 0) return {pref, "Error parsing request: count must be greter then zero"};
if(pdims[ind].beg >= dlen) return {pref, MString("Error parsing request: start index ") + pdims[ind].beg + " must be lesser then " + pdims[ind].name + " size " + dlen};
if(pdims[ind].beg + pdims[ind].count > dlen)
return {pref, MString("Error parsing request: start index ") + pdims[ind].beg + " with count " + pdims[ind].count + " exceeds " + pdims[ind].name + " size " + dlen};
// Ignore hyperplanes in requests for calculation of data dimensionality
if(pdims[transindex.back()].count == 1) transindex.pop_back();
}
if(transindex.size() != Dimensionity<Data>())
return {pref, MString("Output data dimensions (") + Dimensionity<Data>() + ") not corresponding request dimensions (" + transindex.size() + ")"};
switch(C::VarT(vname))
{
case(C::VarType::UNDEF): return {pref, "No variable with name " + vname + " (impossible)"};
case(C::VarType::FLOAT): return Read<typename C::template Type<C::VarType::FLOAT>>(vname, transindex, data, transform, pdims);
case(C::VarType::DOUBLE): return Read<typename C::template Type<C::VarType::DOUBLE>>(vname, transindex, data, transform, pdims);
case(C::VarType::INT1): return Read<typename C::template Type<C::VarType::INT1>>(vname, transindex, data, transform, pdims);
case(C::VarType::INT2): return Read<typename C::template Type<C::VarType::INT2>>(vname, transindex, data, transform, pdims);
case(C::VarType::INT4): return Read<typename C::template Type<C::VarType::INT4>>(vname, transindex, data, transform, pdims);
case(C::VarType::INT8): return Read<typename C::template Type<C::VarType::INT8>>(vname, transindex, data, transform, pdims);
case(C::VarType::UINT1): return Read<typename C::template Type<C::VarType::INT1>>(vname, transindex, data, transform, pdims);
}
return {pref, "Internal error (impossible)"};
}
// Request by string argument
template<class Data, class Transform> Error Read(const MString& vname, Data& data, Transform transform, const MString& request) const
{
static const MString pref = "NcZarrRead::Read";
std::vector<struct DimReq> pdims;
// Parse request
const auto dimdesc = request.Split(";, \t");
if(dimdesc.size() == 0) return {pref, "Empty request"};
for(const auto& dd: dimdesc)
{
const auto dimpar = dd.Split(":", true);
if(dimpar.size() == 1) // Only name, so, we request full length
pdims.emplace_back(dimpar[0]);
else if(dimpar.size() == 2) // Name and first index
pdims.emplace_back(dimpar[0], dimpar[1].ToInteger<size_t>());
else if(dimpar.size() == 3) // Name, first index, count
pdims.emplace_back(dimpar[0], dimpar[1].ToInteger<size_t>(), dimpar[2].ToInteger<size_t>());
else
return {pref, "Can't parse expression " + dd};
}
return Read(vname, data, transform, pdims);
}
// Request full one-dimensional variable
template<class Type> Error Read(const MString& vname, std::vector<Type>& out) const
{
const auto& dnames = C::DimNames(vname);
if(dnames.size() > 0) out.resize(C::DimSize(dnames[0]));
auto data = [&vec = out](size_t i) -> Type& { return vec[i]; };
return Read(vname, data, std::identity());
}
};

175
include/zarr.h

@ -0,0 +1,175 @@
#pragma once
#include "GPL.h"
#include "cache.h"
#include "curlfuncs.h"
#include "nczarrcommon.h"
#include <json/json.h>
#include <variant>
class ZarrTypes: public NcZarrTypes
{
protected:
template<class VType> class ReadedData
{
//public:
using Vec = std::vector<size_t>;
private:
Vec start, chunkstart;
ArrCounter mainind, chunkind, inchunkind;
std::vector<std::unique_ptr<VType[]>> data;
public:
ReadedData(): mainind(Vec()), chunkind(Vec()), inchunkind(Vec()) {}
ReadedData(size_t N, const size_t* start, const size_t* count, const size_t* csize, std::vector<std::unique_ptr<VType[]>>&& d):
start(start, start + N),
chunkstart(
[](size_t N, const size_t* st, const size_t* cs)
{
Vec out(N);
for(size_t i = 0; i < N; i++) out[i] = st[i] / cs[i];
return out;
}(N, start, csize)),
mainind(Vec(count, count + N)),
chunkind(
[](size_t N, const size_t* st, const size_t* cn, const size_t* cs)
{
Vec out(N);
for(size_t i = 0; i < N; i++) out[i] = (st[i] + cn[i]) / cs[i] - st[i] / cs[i] + 1;
return out;
}(N, start, count, csize)),
inchunkind(Vec(csize, csize + N)),
data(std::move(d))
{
}
VType operator()(size_t lini) const
{
Vec ind = mainind.Index(lini, mainind.Count());
Vec cind(ind.size()), inind(ind.size());
for(size_t i = 0; i < ind.size(); i++)
{
cind[i] = (ind[i] + start[i]) / inchunkind.Count(i) - chunkstart[i]; // indes of chunk
inind[i] = (ind[i] + start[i]) % inchunkind.Count(i); // index inside chunk
}
size_t chunk = chunkind.Index(cind);
size_t inside = inchunkind.Index(inind);
return data[chunk][inside];
}
};
private:
// Create attribute from json value
static AttVT CreateAtt(const Json::Value& val)
{
if(val.type() == Json::intValue) return AttVT{std::in_place_type<int8>, val.asInt64()};
if(val.type() == Json::uintValue) return AttVT(std::in_place_type<uint8>, val.asUInt64());
if(val.type() == Json::realValue) return AttVT(std::in_place_type<double>, val.asDouble());
if(val.type() == Json::stringValue)
{
auto str = val.asString();
return AttVT(std::in_place_type<MString>, MString(str.c_str(), str.size()));
}
if(val.type() == Json::booleanValue) return AttVT(std::in_place_type<bool>, val.asBool());
return AttVT();
}
public:
// Read attributes from .zattrs
static auto ReadAtts(const Json::Value& obj)
{
std::vector<Attribute> out;
if(obj.type() != Json::objectValue) return out;
const auto keys = obj.getMemberNames();
for(const auto& key: keys)
if(key != "_ARRAY_DIMENSIONS") out.emplace_back(key, CreateAtt(obj[key]));
return out;
}
};
class ZarrFunctions: public ZarrTypes
{
std::unique_ptr<GenericCache> cache;
CURLRAII chandle;
MString url;
std::vector<std::vector<size_t>> chunks;
// Find variable names in metadata
static std::vector<MString> ReadVarNames(const Json::Value& meta);
Error AddVar(const MString& name, const Json::Value& zattrs, const Json::Value& zarray);
protected:
ZarrFunctions()
{
auto oldprefix = michlib::GPL.UsePrefix("ZARR");
cache.reset(CreateCache(michlib::GPL.ParameterSValue("Cache", "")));
michlib::GPL.UsePrefix(oldprefix);
if(!cache)
{
michlib::errmessage("Can't init data cache");
cache.reset(new FakeCache);
}
}
template<class VType> RetVal<ReadedData<VType>> Read(const MString& var, const size_t* start, const size_t* count) const
{
using Vec = std::vector<size_t>;
size_t ind = FindInd(var, vars);
const size_t N = vars[ind].NDim();
const auto& csize = chunks[ind];
Vec chunkstart(
[](size_t N, const size_t* st, const size_t* cs)
{
Vec out(N);
for(size_t i = 0; i < N; i++) out[i] = st[i] / cs[i];
return out;
}(N, start, csize.data()));
ArrCounter chunkind(
[](size_t N, const size_t* st, const size_t* cn, const size_t* cs)
{
Vec out(N);
for(size_t i = 0; i < N; i++) out[i] = (st[i] + cn[i] - 1) / cs[i] - st[i] / cs[i] + 1;
return out;
}(N, start, count, csize.data()));
bool havefill = vars[ind].Fill().index() > 0;
VType fill = std::visit(
[](auto v)
{
if constexpr(std::is_convertible_v<decltype(v), VType>)
return static_cast<VType>(v);
else
return std::numeric_limits<VType>::max();
},
vars[ind].Fill());
std::vector<std::unique_ptr<VType[]>> cdata;
size_t chunksize = 1;
for(const auto c: csize) chunksize *= c;
cdata.resize(chunkind.N());
for(; chunkind; ++chunkind)
{
cdata[chunkind.Index()].reset(new VType[chunksize]);
auto res = GetChunk(var, chunkind.VIndex(chunkstart), chunksize, sizeof(VType), cdata[chunkind.Index()].get(), havefill ? &fill : nullptr);
if(!res) return res;
}
return ReadedData<VType>(N, start, count, csize.data(), std::move(cdata));
}
Error GetChunk(const MString& var, const std::vector<size_t>& chunkind, size_t chunksize, size_t elsize, void* data, const void* fill) const;
public:
Error Open(const MString& product, const MString& dataset, bool time = true);
};
using Zarr = NcZarrRead<ZarrFunctions>;

2
michlib

@ -1 +1 @@
Subproject commit e2882902b88229bb4b0e6fbeb76c79ac6d46d53d
Subproject commit f380988909fadd7f42c7ce09288c4ff3198ca318

6
sources/AVISO.h

@ -1,7 +1,7 @@
#pragma once
#include "layereddata.h"
#include "layereddataz.h"
class AVISOData: public LayeredData
class AVISOData: public LayeredDataZ
{
enum Type
{
@ -49,6 +49,6 @@ class AVISOData: public LayeredData
return "Unknown dataset: " + dataset;
SetTitle(DataTitle());
return LayeredData::Open(dataset);
return LayeredDataZ::Open(dataset);
}
};

2
sources/AVISOLOCAL.h

@ -49,7 +49,7 @@ class AVISOLOCALData: public NCFuncs
return times[i];
}
time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
time_t Timestep() const { return isOk() ? (times[1] - times[0]).Seconds() : 0; }
explicit operator bool() const { return times.size() > 0; }

2
sources/BINFILE.h

@ -38,7 +38,7 @@ class BINFILEData
return times[i];
}
time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
time_t Timestep() const { return isOk() ? (times[1] - times[0]).Seconds() : 0; }
explicit operator bool() const { return times.size() > 0; }

146
sources/COPERNICUS.cpp

@ -6,50 +6,9 @@
using michlib::GPL;
const MString COPERNICUSData::caturl = "https://stac.marine.copernicus.eu/metadata/catalog.stac.json";
std::pair<Json::Value, MString> COPERNICUSData::GetJSON(const MString& url)
{
Json::Reader reader;
Json::Value obj;
MString content;
auto [val, suc] = cache->Get(url);
if(suc)
content = std::move(val);
else
{
michlib::message(url + " not found in cache, downloading");
auto [out, res] = GetUrl(chandle, url);
if(res != CURLE_OK) return {obj, MString("Can't download JSON: ") + curlerr};
cache->Put(url, out, 3600);
content = std::move(out);
}
reader.parse(content.Buf(), content.Buf() + content.Len(), obj, false);
return {obj, ""};
}
MString COPERNICUSData::ReadURL(const Json::Value& cat, const MString& prod)
{
const auto& links = cat["links"];
if(links.type() != Json::arrayValue) return "";
for(Json::ArrayIndex i = 0; i < links.size(); i++)
{
const auto& titl = links[i]["title"];
const auto& href = links[i]["href"];
if(titl.type() == Json::stringValue && href.type() == Json::stringValue)
{
MString str(titl.asString().c_str());
if(str == prod) return MString(href.asString().c_str());
}
}
return "";
}
std::pair<std::vector<struct FileInfo>, MString> COPERNICUSData::ReadRemoteFileList(const MString& url)
RetVal<std::vector<struct FileInfo>> COPERNICUSData::ReadRemoteFileList(const MString& url) const
{
const static MString pref = "COPERNICUSData::ReadRemoteFileList";
LIBXML_TEST_VERSION
std::vector<struct FileInfo> out;
@ -68,7 +27,7 @@ std::pair<std::vector<struct FileInfo>, MString> COPERNICUSData::ReadRemoteFileL
break;
}
}
if(pos == url.Len()) return {out, "Can't parse url: " + url};
if(pos == url.Len()) return {pref, "Can't parse url: " + url};
bucket = url.SubStr(1, pos);
prefix = url.SubStr(pos + 2, url.Len() - pos - 1);
@ -77,6 +36,7 @@ std::pair<std::vector<struct FileInfo>, MString> COPERNICUSData::ReadRemoteFileL
MString cont;
bool next = true;
CURLRAII chandle;
while(next)
{
MString url = bucket + "?list-type=2&prefix=" + prefix;
@ -84,20 +44,20 @@ std::pair<std::vector<struct FileInfo>, MString> COPERNICUSData::ReadRemoteFileL
cont = "";
auto [data, res] = GetUrl(chandle, url);
if(res != CURLE_OK) return {out, MString("Can't download ") + url + ": " + curlerr};
if(res != CURLE_OK) return {pref, MString("Can't download ") + url + ": " + chandle.Err()};
xmlDocPtr doc = xmlReadMemory(data.Buf(), data.Len(), "data.xml", nullptr, 0);
if(doc == nullptr) return {out, MString("Can't download ") + url + ": XML parse error"};
if(doc == nullptr) return {pref, MString("Can't download ") + url + ": XML parse error"};
auto cur = xmlDocGetRootElement(doc);
if(cur == nullptr)
{
xmlFreeDoc(doc);
return {out, MString("Can't download ") + url + ": empty XML"};
return {pref, MString("Can't download ") + url + ": empty XML"};
}
if(xmlStrEqual(cur->name, (const xmlChar*)"ListBucketResult") == 0)
{
xmlFreeDoc(doc);
return {out, MString("Can't download ") + url + ": unknown XML"};
return {pref, MString("Can't download ") + url + ": unknown XML"};
}
for(const auto* n = cur->children; n; n = n->next)
@ -142,87 +102,49 @@ std::pair<std::vector<struct FileInfo>, MString> COPERNICUSData::ReadRemoteFileL
}
std::sort(out.begin(), out.end(), [](const struct FileInfo& a, const struct FileInfo& b) { return a.name < b.name; });
return {out, ""};
return out;
}
MString COPERNICUSData::Mirror(const CLArgs& args)
Error COPERNICUSData::Mirror(const CLArgs& args) const
{
const static MString pref = "COPERNICUSData::Mirror";
GPL.UsePrefix("COPERNICUS");
// Local directory
MString mirrorroot = GPL.ParameterSValue("MirrorTo", "");
if(!mirrorroot.Exist()) return "Local mirror directory not specified";
// Cache
cache.reset(CreateCache(GPL.ParameterSValue("Cache", "")));
if(!cache)
{
michlib::errmessage("Can't init cache");
cache.reset(new FakeCache);
}
curl_easy_setopt(chandle, CURLOPT_ERRORBUFFER, curlerr);
if(!args.contains("product")) return "Copernicus product not specified";
MString prod = args.at("product");
Json::Value product;
MString produrl;
// Get catalog
{
auto [cat, err] = GetJSON(caturl);
if(err.Exist()) return "Can't download catalog: " + err;
if(cat["title"].type() != Json::stringValue || cat["title"].asString() != "Copernicus Marine Data Store") return "Can't parse catalog";
catalog = std::move(cat);
}
if(!mirrorroot.Exist()) return {pref, "Local mirror directory not specified"};
// Get product
{
auto url = ReadURL(catalog, prod);
if(!url.Exist()) return "Url for product " + prod + " not found in catalog";
produrl = DirName(caturl) + "/" + url;
auto [pr, err] = GetJSON(produrl);
if(err.Exist()) return "Can't download product information from " + produrl + ": " + err;
product = std::move(pr);
}
if(!args.contains("product")) return {pref, "Copernicus product not specified"};
MString prod = args.at("product");
CopernicusCatalog cat;
std::vector<MString> dsets;
if(args.contains("dataset"))
dsets.push_back(args.at("dataset"));
else
{
const auto& links = product["links"];
if(links.type() != Json::arrayValue) return "Can't find information about datasets";
for(Json::ArrayIndex i = 0; i < links.size(); i++)
{
const auto& rel = links[i]["rel"];
const auto& titl = links[i]["title"];
if(rel.type() == Json::stringValue && titl.type() == Json::stringValue && rel.asString() == "item") dsets.push_back(titl.asString().c_str());
}
auto dlist = cat.DatasetList(prod);
if(!dlist) return dlist.Add(pref, "Can't get list of datasets");
dsets = dlist.Value();
}
CURLRAII dhandle;
CURLRAII chandle;
for(const auto& dset: dsets)
{
michlib::message("Mirroring " + dset);
auto url = ReadURL(product, dset);
if(!url.Exist()) return "Url for dataset " + dset + " not found in product description";
MString dseturl = DirName(produrl) + "/" + url;
auto [ds, err] = GetJSON(dseturl);
if(err.Exist()) return "Can't download dataset information from " + dseturl + ": " + err;
const auto& href = ds["assets"]["native"]["href"];
if(href.type() != Json::stringValue) return "Can't find data for dataset " + dset + " from product " + prod;
url = href.asString().c_str();
auto url = cat.DatasetNativeURL(prod, dset);
if(!url) return {pref, "Can't find data for dataset " + dset + " from product " + prod};
MString locroot = mirrorroot + "/" + prod + "/" + dset;
auto [lfiles, lerr] = ReadLocalFileList(locroot);
if(lerr.Exist()) return lerr;
auto lfilesret = ReadLocalFileList(locroot);
if(!lfilesret) return lfilesret.Add(pref, "Can't get local file list");
const auto& lfiles = lfilesret.Value();
auto [rfiles, rerr] = ReadRemoteFileList(url);
if(rerr.Exist()) return rerr;
auto rfilesret = ReadRemoteFileList(url.Value());
if(!rfilesret) return rfilesret.Add(pref, "Can't get remote file list");
const auto& rfiles = rfilesret.Value();
std::vector<size_t> down, rem;
std::vector<std::pair<size_t, size_t>> upd;
@ -259,25 +181,25 @@ MString COPERNICUSData::Mirror(const CLArgs& args)
for(size_t i = 0; i < down.size(); i++)
{
size_t ri = down[i];
auto err = DownloadFile(dhandle, rfiles[ri], locroot);
if(err.Exist()) return err;
auto err = DownloadFile(chandle, rfiles[ri], locroot);
if(!err) return err.Add(pref, "Can't download file");
}
for(size_t i = 0; i < rem.size(); i++)
{
size_t li = rem[i];
auto err = RemoveFile(lfiles[li]);
if(err.Exist()) return err;
if(!err) return err.Add(pref, "Can't remove file");
}
for(size_t i = 0; i < upd.size(); i++)
{
size_t ri = upd[i].first;
size_t li = upd[i].second;
auto err = UpdateFile(dhandle, rfiles[ri], lfiles[li], locroot);
if(err.Exist()) return err;
auto err = UpdateFile(chandle, rfiles[ri], lfiles[li], locroot);
if(!err) return err.Add(pref, "Can't update file");
}
}
return "";
return Error();
}

21
sources/COPERNICUS.h

@ -1,30 +1,15 @@
#pragma once
#include "ParseArgs.h"
#include "cache.h"
#include "curlfuncs.h"
#include "copcat.h"
#include "mdatetime.h"
#include <json/json.h>
using michlib::MDateTime;
using michlib::MString;
class COPERNICUSData
{
static const MString caturl;
std::unique_ptr<GenericCache> cache;
CURLRAII chandle;
Json::Value catalog;
char curlerr[CURL_ERROR_SIZE];
// Get url for product or dataset from catalog
static MString ReadURL(const Json::Value& cat, const MString& prod);
// Download JSON from url
std::pair<Json::Value, MString> GetJSON(const MString& url);
// Get remote file list from url
std::pair<std::vector<struct FileInfo>,MString> ReadRemoteFileList(const MString& url);
RetVal<std::vector<struct FileInfo>> ReadRemoteFileList(const MString& url) const;
public:
static constexpr const char* name = "COPERNICUS";
@ -32,5 +17,5 @@ class COPERNICUSData
COPERNICUSData() = default;
// Main mirror function
MString Mirror(const CLArgs& args);
Error Mirror(const CLArgs& args) const;
};

20
sources/NEMO.h

@ -1,7 +1,7 @@
#pragma once
#include "layereddata.h"
#include "layereddataz.h"
class NEMOData: public LayeredData
class NEMOData: public LayeredDataZ
{
enum Type
{
@ -16,13 +16,10 @@ class NEMOData: public LayeredData
TYPE_BLKSEANRT,
TYPE_MEDSEADT,
TYPE_MEDSEANRT,
TYPE_MEDSEANRT1,
TYPE_BISCDT,
TYPE_BISCNRT,
TYPE_ENWSDT,
TYPE_ENWSNRT,
TYPE_BRITNRT,
TYPE_BRITNRT1
TYPE_ENWSNRT
};
Type type = TYPE_UNKNOWN;
@ -46,13 +43,10 @@ class NEMOData: public LayeredData
case(TYPE_BLKSEANRT): return "NEMO Near-real time time, Black Sea region, daily mean (BLKSEANRT)";
case(TYPE_MEDSEADT): return "NEMO Delayed time, Mediterranean Sea region, daily mean (MEDSEADT)";
case(TYPE_MEDSEANRT): return "NEMO Near-real time, Mediterranean Sea region, daily mean (MEDSEANRT)";
case(TYPE_MEDSEANRT1): return "NEMO Near-real time, Mediterranean Sea region, hourly mean (MEDSEANRT1)";
case(TYPE_BISCDT): return "NEMO Delayed time, Atlantic-Iberian Biscay Irish region, daily mean (BISCDT)";
case(TYPE_BISCNRT): return "NEMO Near-real time, Atlantic-Iberian Biscay Irish region, daily mean (BISCNRT)";
case(TYPE_ENWSDT): return "NEMO Delayed time, Atlantic - European North West Shelf region, daily mean (ENWSDT)";
case(TYPE_ENWSNRT): return "NEMO Near-real time, Atlantic - European North West Shelf region, daily mean (ENWSNRT)";
case(TYPE_BRITNRT): return "NEMO Near-real time, British Islands region, daily mean (BRITNRT)";
case(TYPE_BRITNRT1): return "NEMO Near-real time, British Islands region, 1h resolution (BRITNRT1)";
default: return "No title";
}
}
@ -83,8 +77,6 @@ class NEMOData: public LayeredData
type = TYPE_MEDSEADT;
else if(dataset == "MEDSEANRT")
type = TYPE_MEDSEANRT;
else if(dataset == "MEDSEANRT1")
type = TYPE_MEDSEANRT1;
else if(dataset == "BISCDT")
type = TYPE_BISCDT;
else if(dataset == "BISCNRT")
@ -93,14 +85,10 @@ class NEMOData: public LayeredData
type = TYPE_ENWSDT;
else if(dataset == "ENWSNRT")
type = TYPE_ENWSNRT;
else if(dataset == "BRITNRT")
type = TYPE_BRITNRT;
else if(dataset == "BRITNRT1")
type = TYPE_BRITNRT1;
else
return "Unknown dataset: " + dataset;
SetTitle(DataTitle());
return LayeredData::Open(dataset);
return LayeredDataZ::Open(dataset);
}
};

6
sources/NEMOBIO.h

@ -1,7 +1,7 @@
#pragma once
#include "layereddata.h"
#include "layereddataz.h"
class NEMOBIOData: public LayeredData
class NEMOBIOData: public LayeredDataZ
{
enum Type
{
@ -43,6 +43,6 @@ class NEMOBIOData: public LayeredData
return "Unknown dataset: " + dataset;
SetTitle(DataTitle());
return LayeredData::Open(dataset);
return LayeredDataZ::Open(dataset);
}
};

10
sources/VYLET.cpp

@ -258,7 +258,7 @@ VYLETData::Data VYLETData::ReadL() const
{
time = R2Time((*vylet)[tcol - 1][iy * lons.size() + ix]);
lambda = (*vylet)[lcol - 1][iy * lons.size() + ix];
days = static_cast<real>(time - start) / MDateTime::secondsperday;
days = (time - start).D();
out(ix, iy) = lambda / days;
}
@ -284,7 +284,7 @@ VYLETData::Data VYLETData::ReadT() const
const real acc = vylet->ParameterRValue("accuracy", 1.0);
const real tstep = 2.0 * M_PI * 1000.0 / acc;
const real maxdays = static_cast<real>(end - start) / MDateTime::secondsperday;
const real maxdays = (end - start).D();
MDateTime time;
real days;
@ -297,7 +297,7 @@ VYLETData::Data VYLETData::ReadT() const
x = (*vylet)[2][iy * lons.size() + ix];
y = (*vylet)[3][iy * lons.size() + ix];
time = R2Time((*vylet)[tcol - 1][iy * lons.size() + ix]);
days = static_cast<real>(time - start) / MDateTime::secondsperday;
days = (time - start).D();
if(days <= tstep * 1.5) days = 0.0;
inside = x > xl && x < xr && y > yd && y < yu;
maxtime = days >= maxdays - tstep * 0.5;
@ -505,7 +505,7 @@ VYLETData::Data VYLETData::ReadTmask() const
const real acc = vylet->ParameterRValue("accuracy", 1.0);
const real tstep = 2.0 * M_PI * 1000.0 / acc;
const real maxdays = static_cast<real>(end - start) / MDateTime::secondsperday;
const real maxdays = (end - start).D();
MDateTime time;
real days;
@ -515,7 +515,7 @@ VYLETData::Data VYLETData::ReadTmask() const
for(size_t ix = 0; ix < lons.size(); ix++)
{
time = R2Time((*vylet)[tcol - 1][iy * lons.size() + ix]);
days = static_cast<real>(time - start) / MDateTime::secondsperday;
days = (time - start).D();
maxtime = days >= maxdays - tstep * 0.5;
out(ix, iy) = maxtime ? 1.0 : NAN;
}

6
src/CMakeLists.txt

@ -9,16 +9,18 @@ find_package(CURL REQUIRED)
find_package(LibXml2 REQUIRED)
find_package(SQLite3 REQUIRED)
pkg_check_modules(JSONCPP REQUIRED jsoncpp)
pkg_check_modules(BLOSC REQUIRED blosc)
pkg_check_modules(LIBPQ REQUIRED libpq)
set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
include_directories(${JSONCPP_INCLUDE_DIRS} ${LIBXML2_INCLUDE_DIRS} ${SQLite3_INCLUDE_DIRS})
include_directories(${JSONCPP_INCLUDE_DIRS} ${BLOSC_INCLUDE_DIRS} ${LIBPQ_INCLUDE_DIRS} ${LIBXML2_INCLUDE_DIRS} ${SQLite3_INCLUDE_DIRS})
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} OpenMP::OpenMP_CXX CURL::libcurl ${JSONCPP_LINK_LIBRARIES} LibXml2::LibXml2 SQLite::SQLite3 teos)
target_link_libraries(${EXENAME} ${linker_options} ${netcdf} 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})

168
src/copcat.cpp

@ -0,0 +1,168 @@
#define MICHLIB_NOSOURCE
#include "copcat.h"
#include "GPL.h"
#include "mirrorfuncs.h"
const MString CopernicusCatalog::caturl = "https://stac.marine.copernicus.eu/metadata/catalog.stac.json";
CopernicusCatalog::CopernicusCatalog()
{
// Cache
auto oldprefix = michlib::GPL.UsePrefix("COPERNICUS");
cache.reset(CreateCache(michlib::GPL.ParameterSValue("Cache", "")));
michlib::GPL.UsePrefix(oldprefix);
if(!cache)
{
michlib::errmessage("Can't init cache");
cache.reset(new FakeCache);
}
GetCatalog();
}
Error CopernicusCatalog::GetCatalog()
{
if(Valid()) return Error();
auto ret = GetJSON(caturl);
if(ret)
catalog = ret.Value();
else
return ret.Add("CopernicusCatalog::GetCatalog", "can't download catalog");
return Error();
}
RetVal<std::vector<MString>> CopernicusCatalog::ProductList() const
{
static const MString pref = "CopernicusCatalog::ProductList";
if(!Valid()) return {pref, "no catalog"};
const auto& links = catalog["links"];
if(links.type() != Json::arrayValue) return {pref, "no \"links\" section in the catalog"};
std::vector<MString> out;
for(Json::ArrayIndex i = 0; i < links.size(); i++)
{
const auto& rel = links[i]["rel"];
const auto& titl = links[i]["title"];
if(rel.type() == Json::stringValue && titl.type() == Json::stringValue && rel.asString() == "child") out.emplace_back(titl.asString().c_str());
}
return out;
}
RetVal<MString> CopernicusCatalog::ProductURL(const MString& prod) const
{
static const MString pref = "CopernicusCatalog::ProductURL";
if(!Valid()) return {pref, "no catalog"};
const auto& links = catalog["links"];
if(links.type() != Json::arrayValue) return {pref, "no \"links\" section in the catalog"};
for(Json::ArrayIndex i = 0; i < links.size(); i++)
{
const auto& titl = links[i]["title"];
const auto& href = links[i]["href"];
if(titl.type() == Json::stringValue && href.type() == Json::stringValue && titl.asString().c_str() == prod) return DirName(caturl) + "/" + MString(href.asString().c_str());
}
return {pref, "unknown product: " + prod};
}
RetVal<std::vector<MString>> CopernicusCatalog::DatasetList(const MString& prod) const
{
static const MString pref = "CopernicusCatalog::DatasetList";
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& links = ret.Value()["links"];
if(links.type() != Json::arrayValue) return {pref, "no \"links\" section in the product " + prod + " description"};
std::vector<MString> out;
for(Json::ArrayIndex i = 0; i < links.size(); i++)
{
const auto& rel = links[i]["rel"];
const auto& titl = links[i]["title"];
if(rel.type() == Json::stringValue && titl.type() == Json::stringValue && rel.asString() == "item") out.emplace_back(titl.asString().c_str());
}
return out;
}
RetVal<MString> CopernicusCatalog::DatasetURL(const MString& prod, const MString& dataset) const
{
static const MString pref = "CopernicusCatalog::DatasetURL";
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& links = ret.Value()["links"];
if(links.type() != Json::arrayValue) return {pref, "no \"links\" section in the product " + prod + " description"};
for(Json::ArrayIndex i = 0; i < links.size(); i++)
{
const auto& titl = links[i]["title"];
const auto& href = links[i]["href"];
if(titl.type() == Json::stringValue && href.type() == Json::stringValue && titl.asString().c_str() == dataset) return DirName(url) + "/" + MString(href.asString().c_str());
}
return {pref, "unknown dataset: " + dataset};
}
RetVal<MString> CopernicusCatalog::AssetURL(const MString& prod, const MString& dataset, const MString& asset) const
{
static const MString pref = "CopernicusCatalog::AssetURL";
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& href = ret.Value()["assets"][asset.Buf()]["href"];
if(!href || href.type() != Json::stringValue) return {pref, "href for the asset " + asset + " not found"};
return MString(href.asString().c_str());
}
RetVal<Json::Value> CopernicusCatalog::GetJSON(const MString& url) const
{
const static MString pref = "CopernicusCatalog::GetJSON";
Json::Reader reader;
Json::Value obj;
MString content;
auto [val, suc] = cache->Get(url);
if(suc)
content = std::move(val);
else
{
michlib::message(url + " not found in cache, downloading");
auto [out, res] = GetUrl(chandle, url);
if(res != CURLE_OK) return Error(pref, MString("can't download JSON: ") + chandle.Err());
cache->Put(url, out, 3600);
content = std::move(out);
}
reader.parse(content.Buf(), content.Buf() + content.Len(), obj, false);
return obj;
}

298
src/layereddataz.cpp

@ -0,0 +1,298 @@
#define MICHLIB_NOSOURCE
#include "layereddataz.h"
MString LayeredDataZ::Info() const
{
if(!isOk()) return "";
MString d;
for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")";
std::set<MString> vars;
for(const auto& f: nc) GetVars(f, vars);
MString svars;
{
bool first = true;
for(const auto& v: vars)
{
svars += (first ? "" : ", ") + v;
first = false;
}
}
// clang-format off
return
"Dataset: " + Title() + "\n" +
" Begin date: " + Time(0).ToString() + "\n" +
" End date: " + Time(NTimes()-1).ToString() + "\n" +
" Time step: " + Timestep() + " seconds\n" +
" Time moments: " + NTimes() + "\n" +
" Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" +
" Grid: " + dname.nx + "x" + dname.ny + " (" + lonstep + " x " + latstep + ")\n" +
" Depths:" + d + "\n" +
" Supported variables: " + svars;
// clang-format on
}
MString LayeredDataZ::Open(const MString& dataset)
{
nc.clear();
MString proxyurl = GPL.ParameterSValue("USEPROXY", "");
if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl);
nc.clear();
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() == 0 || words.size() > 2)
{
nc.clear();
return "Invalid url " + url;
}
MString product = words[0];
MString dataset = words.size() == 2 ? words[1] : "";
nc.emplace_back();
{
auto ret = nc.back().Open(product, dataset);
if(!ret)
{
nc.clear();
return "Can't open " + dataset + " of " + product;
}
}
}
else
break;
i++;
}
if(nc.size() == 0) return "No urls for dataset " + dataset + " specified in config";
dname = GetDNames(nc[0]);
if(!(dname.lonname.Exist() && dname.latname.Exist()))
{
nc.clear();
return "Can't find longitude/latitude";
}
if(!dname.timename.Exist())
{
nc.clear();
return "Can't find time";
}
auto cn = GetCNames(nc[0]);
// Read times
for(auto& f: nc)
{
auto ret = f.ReadTimes(cn.timename);
if(!ret)
{
nc.clear();
return "Can't read times";
}
times.insert(times.end(), f.Times().begin(), f.Times().end());
}
std::sort(times.begin(), times.end());
auto last = std::unique(times.begin(), times.end());
times.erase(last, times.end());
depthinv = false;
if(cn.depthname.Exist())
{
auto ret = nc[0].Read(cn.depthname, depths);
if(!ret)
{
nc.clear();
return "Can't read depths";
}
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]);
}
}
else // Surface only data
{
depths.resize(1);
depths[0] = 0;
}
std::vector<double> lons, lats;
{
auto ret = nc[0].Read(cn.lonname, lons);
if(!ret)
{
nc.clear();
return "Can't get longitudes";
}
}
{
auto ret = nc[0].Read(cn.latname, lats);
if(!ret)
{
nc.clear();
return "Can't get latitudes";
}
}
lonb = lons[0];
latb = lats[0];
lone = lons.back();
late = lats.back();
lonstep = (lone - lonb) / (dname.nx - 1);
latstep = (late - latb) / (dname.ny - 1);
return "";
}
std::pair<const BaseParameters*, MString> LayeredDataZ::Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const
{
std::unique_ptr<struct Parameters> ppar{new struct Parameters};
ppar->layer = args.contains("layer") ? args.at("layer").ToInteger<size_t>() : 0;
if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"};
real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer);
{
auto dom = DetGeoDomain(lonb, lone);
real lon1 = ToGeoDomain(reg.lonb, dom);
real lon2 = ToGeoDomain(reg.lone, dom);
real lat1 = reg.latb;
real lat2 = reg.late;
bool global = lone - lonb + 1.5 * lonstep > 360.0;
// Special case when the longitude lies in a small sector between the end and the start
if(global)
{
if(lon1 < lonb) lon1 = lone;
if(lon2 > lone) lon2 = lonb;
}
else
{
if(lon1 < lonb) lon1 = lonb;
if(lon2 > lone) lon2 = lone;
}
ppar->xb = static_cast<size_t>(Floor((lon1 - lonb) / lonstep));
ppar->xe = static_cast<size_t>(Ceil((lon2 - lonb) / lonstep));
if(ppar->xb == ppar->xe) return {nullptr, "Lonb must be not equal late"};
if(!global && ppar->xb > ppar->xe) return {nullptr, "Lonb must be lesser then lone"};
ppar->yb = static_cast<size_t>(Floor((lat1 - latb) / latstep));
ppar->ye = static_cast<size_t>(Ceil((lat2 - latb) / latstep));
if(ppar->ye > dname.ny - 1) ppar->ye = dname.ny - 1;
if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"};
if(depth < 0.0 || depth > depths.back())
ppar->layer = (depth < 0.0) ? 0 : (depths.size() - 1);
else
for(size_t i = 0; i < depths.size() - 1; i++)
{
if(depth >= depths[i] && depth <= depths[i + 1])
{
ppar->layer = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1);
break;
}
}
if(depthinv) ppar->layer = depths.size() - ppar->layer - 1;
}
pars.SetParameter("depth", Depth(ppar->layer));
pars.SetParameter("layer", ppar->layer);
pars.SetParameter("dataset", Title());
pars.SetParameter("lonb", Lon(ppar->xb));
pars.SetParameter("latb", Lat(ppar->yb));
pars.SetParameter("lone", Lon(ppar->xe));
pars.SetParameter("late", Lat(ppar->ye));
return {ppar.release(), ""};
}
bool LayeredDataZ::Read(const MString& vname, std::map<MString, LayeredDataZ::Data>& cache, const BaseParameters* ip, size_t i) const
{
if(cache.contains(vname)) return true;
if(!isOk()) return false;
auto p = dynamic_cast<const struct Parameters*>(ip);
auto [name, id, tid] = VarNameLoc(vname, times[i]);
if(!name.Exist()) // Conversion read
return TransformationRead(this, vname, cache, ip, i);
// Direct read
bool nodepth = false;
Data data;
//auto head = nc[id]->Header();
for(const auto& v: nc[id].Vars())
if(v.Name() == name)
{
if(v.NDim() == 3) nodepth = true;
if(v.Type() == NcZarrTypes::VarType::INT2) data = ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
if(v.Type() == NcZarrTypes::VarType::INT4) data = ReadVarRaw<int>(nc[id], name, tid, nodepth, p);
if(v.Type() == NcZarrTypes::VarType::FLOAT) data = ReadVarRaw<float>(nc[id], name, tid, nodepth, p);
if(v.Type() == NcZarrTypes::VarType::DOUBLE) data = ReadVarRaw<double>(nc[id], name, tid, nodepth, p);
if(data)
{
cache[vname] = std::move(data);
return true;
}
}
return false;
}
template<class DataType> LayeredDataZ::Data LayeredDataZ::ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct LayeredDataZ::Parameters* p) const
{
real unitmul = 1.0;
//DataType fill;
real offset = 0.0, scale = 1.0;
if(f.HasAtt(name, "add_offset")) offset = f.AttReal(name, "add_offset");
if(f.HasAtt(name, "scale_factor")) scale = f.AttReal(name, "scale_factor");
MString unit;
if(f.HasAtt(name, "units")) unit = f.AttString(name, "units");
if(unit == "m s-1" || unit == "m/s")
{
unitmul = 100.0;
unit = "cm/s";
}
Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (dname.nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep, std::move(unit));
auto trans = [scale, offset, unitmul](auto raw) -> DataType { return (raw * scale + offset) * unitmul; };
if(p->xb < p->xe)
{
auto ret = nodepth ? f.Read(name, data, trans, {dname.lonname, p->xb, p->xe - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
: f.Read(name, data, trans, {dname.lonname, p->xb, p->xe - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1},
{dname.depthname, p->layer, 1});
if(!ret) return Data();
}
else
{
{
auto ret = nodepth ? f.Read(name, data, trans, {dname.lonname, p->xb, dname.nx - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
: f.Read(name, data, trans, {dname.lonname, p->xb, dname.nx - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1},
{dname.depthname, p->layer, 1});
if(!ret) return Data();
}
{
size_t shift = dname.nx - p->xb + 1;
auto shifteddata = [&data, shift](size_t ix, size_t iy) -> real& { return data(ix + shift, iy); };
auto ret =
nodepth ? f.Read(name, shifteddata, trans, {dname.lonname, 0, p->xe + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
: f.Read(name, shifteddata, trans, {dname.lonname, 0, p->xe + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, p->layer, 1});
if(!ret) return Data();
}
}
return data;
}

60
src/mirrorfuncs.cpp

@ -27,32 +27,34 @@ bool MakePath(const MString& dname)
return true;
}
std::pair<std::vector<struct FileInfo>, MString> ReadLocalFileList(const MString& dir, const MString& path)
RetVal<std::vector<struct FileInfo>> ReadLocalFileList(const MString& dir, const MString& path)
{
const static MString pref = "ReadLocalFileList";
std::vector<struct FileInfo> out;
DIRRAII dhandle;
MakePath(dir);
dhandle.reset(opendir(dir.Buf()));
if(!dhandle) return {out, "Can't open directory " + path + (path.Exist() ? "/" : "") + dir};
if(!dhandle) return {pref, "Can't open directory " + path + (path.Exist() ? "/" : "") + dir};
int dfd = dirfd(dhandle);
errno = 0;
struct dirent* dent = readdir(dhandle);
if(errno != 0) return {out, "Can't read directory " + path + (path.Exist() ? "/" : "") + dir};
if(errno != 0) return {pref, "Can't read directory " + path + (path.Exist() ? "/" : "") + dir};
struct stat st;
do {
if(dent->d_name[0] != '.')
{
int ret = fstatat(dfd, dent->d_name, &st, AT_SYMLINK_NOFOLLOW);
if(ret != 0) return {out, "Can't stat " + path + "/" + dir + "/" + dent->d_name};
if(ret != 0) return {pref, "Can't stat " + path + "/" + dir + "/" + dent->d_name};
if(S_ISDIR(st.st_mode)) // Directory, recurse
{
auto [list, err] = ReadLocalFileList(dir + "/" + dent->d_name, path + (path.Exist() ? "/" : "") + dent->d_name);
if(err.Exist()) return {out, err};
out.insert(out.end(), list.begin(), list.end());
auto list = ReadLocalFileList(dir + "/" + dent->d_name, path + (path.Exist() ? "/" : "") + dent->d_name);
if(!list) return list;
out.insert(out.end(), list.Value().begin(), list.Value().end());
}
if(S_ISREG(st.st_mode)) // Regular file
{
@ -63,25 +65,25 @@ std::pair<std::vector<struct FileInfo>, MString> ReadLocalFileList(const MString
dent = readdir(dhandle);
} while(dent != nullptr || errno != 0);
if(errno != 0) return {out, "Can't read directory " + path + "/" + dir};
if(errno != 0) return {pref, "Can't read directory " + path + "/" + dir};
std::sort(out.begin(), out.end(), [](const struct FileInfo& a, const struct FileInfo& b) { return a.name < b.name; });
return {out, ""};
return out;
}
MString DownloadFile(const CURLRAII& chandle, const struct FileInfo& rinfo, const MString& root)
Error DownloadFile(const CURLRAII& chandle, const struct FileInfo& rinfo, const MString& root)
{
const static MString pref = "DownloadFile";
message("Downloading " + rinfo.url);
MString dname = DirName(rinfo.name), fname = FileName(rinfo.name);
FD fd;
if(!MakePath(root + "/" + dname)) return "Can't create directory " + root + "/" + dname;
if(!MakePath(root + "/" + dname)) return {pref, "Can't create directory " + root + "/" + dname};
fd.Reset(creat((root + "/" + rinfo.name).Buf(), 0644));
if(!fd) return "Can't create file " + root + "/" + rinfo.name;
if(!fd) return {pref, "Can't create file " + root + "/" + rinfo.name};
char errbuf[CURL_ERROR_SIZE];
int cfd = fd.Get();
curl_easy_setopt(chandle, CURLOPT_ERRORBUFFER, errbuf);
int cfd = fd.Get();
curl_easy_setopt(chandle, CURLOPT_WRITEFUNCTION, Write2File);
curl_easy_setopt(chandle, CURLOPT_WRITEDATA, &cfd);
curl_easy_setopt(chandle, CURLOPT_URL, rinfo.url.Buf());
@ -89,7 +91,7 @@ MString DownloadFile(const CURLRAII& chandle, const struct FileInfo& rinfo, cons
if(res != CURLE_OK)
{
unlink((root + "/" + rinfo.name).Buf());
return MString("Can't download file: ") + errbuf;
return {pref, MString("Can't download file: ") + chandle.Err()};
}
{
@ -101,30 +103,32 @@ MString DownloadFile(const CURLRAII& chandle, const struct FileInfo& rinfo, cons
if(ret != 0)
{
unlink((root + "/" + rinfo.name).Buf());
return "Can't set mtime for file: " + root + "/" + rinfo.name;
return {pref, "Can't set mtime for file: " + root + "/" + rinfo.name};
}
}
return "";
return Error();
}
MString RemoveFile(const struct FileInfo& linfo)
Error RemoveFile(const struct FileInfo& linfo)
{
const static MString pref = "RemoveFile";
message("Remove " + linfo.url);
int ret = unlink(linfo.url.Buf());
if(ret != 0) return "Can't remove file " + linfo.url;
return "";
if(ret != 0) return {pref, "Can't remove file " + linfo.url};
return Error();
}
MString UpdateFile(const CURLRAII& chandle, const struct FileInfo& rinfo, const struct FileInfo& linfo, const MString& root)
Error UpdateFile(const CURLRAII& chandle, const struct FileInfo& rinfo, const struct FileInfo& linfo, const MString& root)
{
MString err;
const static MString pref = "UpdateFile";
message("Update " + linfo.url);
err = RemoveFile(linfo);
if(err.Exist()) return err;
err = DownloadFile(chandle, rinfo, root);
if(err.Exist()) return err;
auto rm = RemoveFile(linfo);
if(!rm) return rm.Add(pref, "Can't remove file");
auto df = DownloadFile(chandle, rinfo, root);
if(!df) return df.Add(pref, "Can't download file");
return "";
return Error();
}

104
src/ncfuncs.cpp

@ -31,6 +31,37 @@ NCFuncs::CoordNames NCFuncs::GetDNames(const NCFileA& nc)
return out;
}
template<class NcZarrFunctions> NCFuncs::CoordNames NCFuncs::GetDNames(const NcZarrRead<NcZarrFunctions>& nc)
{
CoordNames out;
for(const auto& dim: nc.Dims())
{
if(dim.Name() == "lon" || dim.Name() == "longitude")
{
out.lonname = dim.Name();
out.nx = dim.Size();
}
if(dim.Name() == "lat" || dim.Name() == "latitude")
{
out.latname = dim.Name();
out.ny = dim.Size();
}
if(dim.Name() == "depth" || dim.Name() == "elevation")
{
out.depthname = dim.Name();
out.nz = dim.Size();
}
if(dim.Name() == "time")
{
out.timename = dim.Name();
out.nt = dim.Size();
}
}
return out;
}
template NCFuncs::CoordNames NCFuncs::GetDNames<ZarrFunctions>(const NcZarrRead<ZarrFunctions>&);
NCFuncs::CoordNames NCFuncs::GetCNames(const NCFileA& nc)
{
CoordNames out;
@ -64,6 +95,46 @@ NCFuncs::CoordNames NCFuncs::GetCNames(const NCFileA& nc)
return out;
}
template<class NcZarrFunctions> NCFuncs::CoordNames NCFuncs::GetCNames(const NcZarrRead<NcZarrFunctions>& nc)
{
CoordNames out;
for(const auto& v: nc.Vars()) // Try to define coordinates by attribute standard_name or attribute axis
{
auto havestname = nc.HasAtt(v.Name(), "standard_name");
auto haveaxis = nc.HasAtt(v.Name(), "axis");
if(!(havestname || haveaxis)) continue;
auto stname = nc.AttString(v.Name(), "standard_name");
auto axis = nc.AttString(v.Name(), "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.lonname.Exist() && axis == "X") islon = true;
if(!out.latname.Exist() && axis == "Y") islat = true;
if(!out.depthname.Exist() && axis == "Z") isdepth = true;
if(!out.timename.Exist() && axis == "T") istime = true;
if(islon) out.lonname = v.Name();
if(islat) out.latname = v.Name();
if(isdepth) out.depthname = v.Name();
if(istime) out.timename = v.Name();
if(islon) out.nx = v.Dims().size();
if(islat) out.ny = v.Dims().size();
if(isdepth) out.nz = v.Dims().size();
if(istime) out.nt = v.Dims().size();
}
// If time not found just check variable "time"
if(!out.timename.Exist() && nc.HasVar("time")) out.timename = "time";
return out;
}
template NCFuncs::CoordNames NCFuncs::GetCNames<ZarrFunctions>(const NcZarrRead<ZarrFunctions>&);
void NCFuncs::GetVars(const NCFileA& nc, std::set<MString>& vars)
{
auto head = nc.Header();
@ -83,6 +154,26 @@ void NCFuncs::GetVars(const NCFileA& nc, std::set<MString>& vars)
if(vars.contains("ssh")) vars.emplace("vgeo");
}
template<class NcZarrFunctions> void NCFuncs::GetVars(const NcZarrRead<NcZarrFunctions>& nc, std::set<MString>& vars)
{
for(const auto& v: nc.Vars())
{
if(!nc.HasAtt(v.Name(), "standard_name")) continue;
auto ret = nc.AttString(v.Name(), "standard_name");
if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret));
}
if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens");
if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp");
if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp");
if(vars.contains("u") && vars.contains("v")) vars.emplace("U");
if(vars.contains("u") && vars.contains("v")) vars.emplace("U2");
if(vars.contains("ssh")) vars.emplace("ugeo");
if(vars.contains("ssh")) vars.emplace("vgeo");
}
template void NCFuncs::GetVars<ZarrFunctions>(const NcZarrRead<ZarrFunctions>&, std::set<MString>&);
std::tuple<MDateTime, time_t, bool> NCFuncs::Refdate(const MString& refdate)
{
MDateTime out;
@ -182,3 +273,16 @@ bool NCFuncs::HaveVar(const NCFileA& nc, const MString& vname)
}
return false;
}
template<class NcZarrFunctions> bool NCFuncs::HaveVar(const NcZarrRead<NcZarrFunctions>& nc, const MString& vname)
{
for(const auto& v: nc.Vars())
{
if(!nc.HasAtt(v.Name(), "standard_name")) continue;
auto stname = nc.AttString(v.Name(), "standard_name");
if(StName2Name(stname) == vname) return true;
}
return false;
}
template bool NCFuncs::HaveVar<ZarrFunctions>(const NcZarrRead<ZarrFunctions>&, const MString&);

231
src/ncsimple.cpp

@ -0,0 +1,231 @@
#define MICHLIB_NOSOURCE
#include "ncsimple.h"
RetVal<std::vector<NCSimpleFunctions::Attribute>> NCSimpleFunctions::ReadAtts(int vid) const
{
static const MString pref = "NCSimple::ReadAtts";
int natt;
int ret;
if(vid == NC_GLOBAL)
ret = nc_inq_natts(ncid, &natt);
else
ret = nc_inq_var(ncid, vid, nullptr, nullptr, nullptr, nullptr, &natt);
if(ret != NC_NOERR) return Error(pref, MString("Can't inquire number of attributes: ") + nc_strerror(ret));
std::vector<Attribute> out;
char name[NC_MAX_NAME + 1];
for(int aid = 0; aid < natt; aid++)
{
nc_type type;
size_t len;
ret = nc_inq_attname(ncid, vid, aid, name);
if(ret != NC_NOERR) return Error(pref, MString("Can't inquire attribute name: ") + nc_strerror(ret));
ret = nc_inq_atttype(ncid, vid, name, &type);
if(ret != NC_NOERR) return Error(pref, MString("Can't inquire attribute type: ") + nc_strerror(ret));
ret = nc_inq_attlen(ncid, vid, name, &len);
if(ret != NC_NOERR) return Error(pref, MString("Can't inquire attribute length: ") + nc_strerror(ret));
if(type == NC_DOUBLE || type == NC_FLOAT)
{
if(len == 1)
{
double d;
ret = nc_get_att_double(ncid, vid, name, &d);
if(ret != NC_NOERR) return Error(pref, MString("Can't read attribute ") + name + ": " + nc_strerror(ret));
out.emplace_back(MString(name), d);
}
else
{
std::vector<double> dd(len);
ret = nc_get_att_double(ncid, vid, name, dd.data());
if(ret != NC_NOERR) return Error(pref, MString("Can't read attribute ") + name + ": " + nc_strerror(ret));
for(size_t i = 0; i < dd.size(); i++) out.emplace_back(MString(name) + "[" + i + "]", dd[i]);
}
}
else if(type == NC_BYTE || type == NC_SHORT || type == NC_INT || type == NC_INT64)
{
if(len == 1)
{
long long i;
ret = nc_get_att_longlong(ncid, vid, name, &i);
if(ret != NC_NOERR) return Error(pref, MString("Can't read attribute ") + name + ": " + nc_strerror(ret));
out.emplace_back(MString(name), int_cast<int8>(i));
}
else
{
std::vector<long long> ii(len);
ret = nc_get_att_longlong(ncid, vid, name, ii.data());
if(ret != NC_NOERR) return Error(pref, MString("Can't read attribute ") + name + ": " + nc_strerror(ret));
for(size_t i = 0; i < ii.size(); i++) out.emplace_back(MString(name) + "[" + i + "]", int_cast<int8>(ii[i]));
}
}
else if(type == NC_UBYTE || type == NC_USHORT || type == NC_UINT || type == NC_UINT64)
{
if(len == 1)
{
unsigned long long u;
ret = nc_get_att_ulonglong(ncid, vid, name, &u);
if(ret != NC_NOERR) return Error(pref, MString("Can't read attribute ") + name + ": " + nc_strerror(ret));
out.emplace_back(MString(name), int_cast<uint8>(u));
}
else
{
std::vector<unsigned long long> uu(len);
ret = nc_get_att_ulonglong(ncid, vid, name, uu.data());
if(ret != NC_NOERR) return Error(pref, MString("Can't read attribute ") + name + ": " + nc_strerror(ret));
for(size_t i = 0; i < uu.size(); i++) out.emplace_back(MString(name) + "[" + i + "]", int_cast<uint8>(uu[i]));
}
}
else if(type == NC_CHAR)
{
std::vector<char> ss(len + 1, 0);
ret = nc_get_att_text(ncid, vid, name, ss.data());
if(ret != NC_NOERR) return Error(pref, MString("Can't read attribute ") + name + ": " + nc_strerror(ret));
out.emplace_back(MString(name), MString(ss.data()));
}
else
return Error(pref, MString("Unsupported type of attribute ") + name);
// Ignore all other types
}
return out;
}
Error NCSimpleFunctions::Open(const MString& filename)
{
static const MString pref = "NCSimple::Open";
// Cleanup
gats.clear();
dims.clear();
vars.clear();
nc_close(ncid);
std::vector<Attribute> newgats;
std::vector<Dimension> newdims;
std::vector<Variable> newvars;
char name[NC_MAX_NAME + 1];
// Open
{
auto ret = nc_open(filename.Buf(), 0, &ncid);
if(ret != NC_NOERR) return Error(pref, "Can't open file " + filename + ": " + nc_strerror(ret));
}
// Dimensions
{
int ndim;
auto ret = nc_inq_dimids(ncid, &ndim, nullptr, 1);
if(ret != NC_NOERR) return Error(pref, "Can't inquire number of dimensions in file " + filename + ": " + nc_strerror(ret));
std::vector<int> dimids(ndim);
ret = nc_inq_dimids(ncid, nullptr, dimids.data(), 1);
if(ret != NC_NOERR) return Error(pref, "Can't inquire dimension ids in file " + filename + ": " + nc_strerror(ret));
size_t len;
for(const auto id: dimids)
{
ret = nc_inq_dim(ncid, id, name, &len);
if(ret != NC_NOERR) return Error(pref, "Can't inquire dimension name and size in file " + filename + ": " + nc_strerror(ret));
newdims.emplace_back(name, len);
}
}
// Global attributes
{
auto ret = ReadAtts(NC_GLOBAL);
if(!ret) return ret.Add(pref, "Can't read global attributes in file " + filename);
newgats = std::move(ret.Value());
}
// Variables
{
int nvar;
auto ret = nc_inq_varids(ncid, &nvar, nullptr);
if(ret != NC_NOERR) return Error(pref, "Can't inquire number of variables in file " + filename + ": " + nc_strerror(ret));
std::vector<int> varids(nvar);
ret = nc_inq_varids(ncid, nullptr, varids.data());
if(ret != NC_NOERR) return Error(pref, "Can't inquire variables ids in file " + filename + ": " + nc_strerror(ret));
for(const auto vid: varids)
{
nc_type nctype;
int ndim;
auto ret = nc_inq_var(ncid, vid, name, &nctype, &ndim, nullptr, nullptr);
if(ret != NC_NOERR) return Error(pref, "Can't inquire variable info in file " + filename + ": " + nc_strerror(ret));
VarType vt = VarType::UNDEF;
if(nctype == NC_FLOAT) vt = VarType::FLOAT;
if(nctype == NC_DOUBLE) vt = VarType::DOUBLE;
if(nctype == NC_BYTE) vt = VarType::INT1;
if(nctype == NC_SHORT) vt = VarType::INT2;
if(nctype == NC_INT) vt = VarType::INT4;
if(nctype == NC_INT64) vt = VarType::INT8;
if(nctype == NC_UBYTE) vt = VarType::UINT1;
if(vt == VarType::UNDEF) return Error(pref, "Unsupported type of variable " + MString(name) + " in file " + filename);
std::vector<int> dimids(ndim);
ret = nc_inq_vardimid(ncid, vid, dimids.data());
if(ret != NC_NOERR) return Error(pref, "Can't inquire variable dimensions in file " + filename + ": " + nc_strerror(ret));
std::vector<size_t> dims;
char dname[NC_MAX_NAME + 1];
for(const auto did: dimids)
{
auto ret = nc_inq_dimname(ncid, did, dname);
if(ret != NC_NOERR) return Error(pref, "Can't inquire dimension name in file " + filename + ": " + nc_strerror(ret));
size_t ind = newdims.size();
for(size_t i = 0; i < newdims.size() && ind == newdims.size(); i++)
if(dname == newdims[i].Name()) ind = i;
if(ind == newdims.size()) return Error(pref, "Can't find dimension " + MString(dname) + " of variable " + name + " in file " + filename);
dims.push_back(ind);
}
auto atts = ReadAtts(vid);
const char fname[] = "_FillValue";
if(!atts) return Error(pref, "Can't get attributes of variable " + MString(name) + " in file " + filename);
std::vector<Attribute> vatts = std::move(atts.Value());
Variable::FillType fill;
if(FindInd(fname, vatts) != vatts.size())
{
if(nctype == NC_FLOAT || nctype == NC_DOUBLE)
{
double d;
auto ret = nc_get_att_double(ncid, vid, fname, &d);
if(ret != NC_NOERR) return Error(pref, "Can't get fill value for the variable " + MString(name) + " in file " + filename + ": " + nc_strerror(ret));
fill = d;
}
if(nctype == NC_BYTE || nctype == NC_SHORT || nctype == NC_INT || nctype == NC_INT64)
{
long long l;
auto ret = nc_get_att_longlong(ncid, vid, fname, &l);
if(ret != NC_NOERR) return Error(pref, "Can't get fill value for the variable " + MString(name) + " in file " + filename + ": " + nc_strerror(ret));
fill = int_cast<int8>(l);
}
if(nctype == NC_UBYTE || nctype == NC_USHORT || nctype == NC_UINT || nctype == NC_UINT64)
{
unsigned long long u;
auto ret = nc_get_att_ulonglong(ncid, vid, fname, &u);
if(ret != NC_NOERR) return Error(pref, "Can't get fill value for the variable " + MString(name) + " in file " + filename + ": " + nc_strerror(ret));
fill = int_cast<uint8>(u);
}
}
newvars.emplace_back(name, vt, std::move(dims), std::move(vatts), fill);
}
}
gats = std::move(newgats);
dims = std::move(newdims);
vars = std::move(newvars);
return Error();
}

249
src/zarr.cpp

@ -0,0 +1,249 @@
#define MICHLIB_NOSOURCE
#include "zarr.h"
#include "copcat.h"
#include <blosc.h>
std::vector<MString> ZarrFunctions::ReadVarNames(const Json::Value& meta)
{
std::vector<MString> out;
if(meta.type() != Json::objectValue) return out;
const auto keys = meta.getMemberNames();
for(const auto& key: keys)
{
if(!key.ends_with("/.zarray")) continue;
const auto vname = key.substr(0, key.size() - 8);
const auto& zattr = meta[vname + "/.zattrs"];
if(!(zattr && zattr.type() == Json::objectValue)) continue;
MString name(vname.c_str(), vname.size());
bool found = false;
for(size_t id = 0; id < out.size(); id++)
if(out[id] == name)
{
found = true;
break;
}
if(!found) out.emplace_back(std::move(name));
}
return out;
}
Error ZarrFunctions::AddVar(const MString& name, const Json::Value& zattrs, const Json::Value& zarray)
{
static const MString pref = "Zarr::AddVar";
VarType newtype;
Variable::FillType fill;
// Checks for parameters in zarray
{
const auto& cid = zarray["compressor"]["id"];
if(!cid || cid.type() != Json::stringValue || cid.asString() != "blosc") return {pref, "Unsupported compressor: " + MString(cid.asString().c_str())};
}
{
const auto& zf = zarray["zarr_format"];
if(!zf || (zf.type() != Json::uintValue && zf.type() != Json::intValue) || zf.asUInt() != 2) return {pref, "Unsupported format version: " + MString(zf.asUInt())};
}
{
const auto& ord = zarray["order"];
if(!ord || ord.type() != Json::stringValue || ord.asString() != "C") return {pref, "Order in not C"};
}
{
const auto& f = zarray["filters"];
if(f.type() != Json::nullValue) return {pref, "Filters is not null"};
}
// Read dtype
{
const auto& dtype = zarray["dtype"];
if(!dtype || dtype.type() != Json::stringValue) return {pref, "No datatype"};
const auto str = dtype.asString();
if(str == "<f4")
newtype = VarType::FLOAT;
else if(str == "<f8")
newtype = VarType::DOUBLE;
else if(str == "|i1")
newtype = VarType::INT1;
else if(str == "|u1")
newtype = VarType::UINT1;
else if(str == "<i2")
newtype = VarType::INT2;
else if(str == "<i4")
newtype = VarType::INT4;
else if(str == "<i8")
newtype = VarType::INT8;
else
return {pref, "Unsupported datatype: " + MString(str.c_str())};
}
// Read fill_value
{
const auto& fillval = zarray["fill_value"];
if(!fillval) // return {pref, "No fillval"};
fill = 0;
else if(fillval.type() == Json::uintValue)
fill = fillval.asUInt64();
else if(fillval.type() == Json::intValue)
fill = fillval.asInt64();
else if(fillval.type() == Json::realValue)
fill = fillval.asDouble();
else if(fillval.type() == Json::stringValue && fillval.asString() == "NaN")
fill = NAN;
}
// Read attributes
auto atts = ReadAtts(zattrs);
std::vector<MString> dnames;
std::vector<size_t> dsizes;
std::vector<size_t> csizes;
std::vector<size_t> dids;
// Read dimensions names
{
const auto& arrdim = zattrs["_ARRAY_DIMENSIONS"];
if(!(arrdim && arrdim.type() == Json::arrayValue)) return {pref, "_ARRAY_DIMENSIONS not found"};
for(Json::ArrayIndex i = 0; i < arrdim.size(); i++)
if(const auto& dim = arrdim[i]; dim.type() == Json::stringValue)
{
const auto val = dim.asString();
dnames.emplace_back(val.c_str(), val.size());
}
}
// Read dimensions sizes
{
const auto& shape = zarray["shape"];
if(!(shape && shape.type() == Json::arrayValue)) return {pref, "shape not found"};
for(Json::ArrayIndex i = 0; i < shape.size(); i++)
if(const auto& s = shape[i]; s.type() == Json::uintValue || s.type() == Json::intValue) dsizes.push_back(s.asUInt());
}
// Read chunk sizes
{
const auto& chunk = zarray["chunks"];
if(!(chunk && chunk.type() == Json::arrayValue)) return {pref, "chunks not found"};
for(Json::ArrayIndex i = 0; i < chunk.size(); i++)
if(const auto& c = chunk[i]; c.type() == Json::uintValue || c.type() == Json::intValue) csizes.push_back(c.asUInt());
}
if(dnames.size() != dsizes.size() || dnames.size() != csizes.size()) return {pref, "shape and chunks are in contradiction"};
dids.resize(dnames.size());
// Check dimensions names and sizes
for(size_t i = 0; i < dnames.size(); i++)
{
bool found = false;
for(size_t id = 0; id < dims.size(); id++)
if(dims[id].Name() == dnames[i])
{
found = true;
if(dims[id].Size() != dsizes[i])
return {pref, "According to previous data, the dimension " + dnames[i] + " has a size of " + dims[id].Size() + ", but here it is defined as " + dsizes[i]};
dids[i] = id;
break;
}
if(!found)
{
dids[i] = dims.size();
dims.emplace_back(dnames[i], dsizes[i]);
}
}
vars.emplace_back(name, newtype, std::move(dids), std::move(atts), fill);
chunks.push_back(std::move(csizes));
return Error();
}
Error ZarrFunctions::GetChunk(const MString& var, const std::vector<size_t>& chunkind, size_t chunksize, size_t elsize, void* data, const void* fill) const
{
static const MString pref = "Zarr::GetChunk";
MString str = url + "/" + var + "/";
for(size_t i = 0; i < chunkind.size(); i++) str += (i == 0 ? "" : ".") + MString(chunkind[i]);
auto [content, suc] = cache->Get(str);
if(!suc)
{
michlib::message(str + " not found in cache, downloading");
CURLRAII myhandle; // TODO: remove this workaround of unknown bug
//auto [out, res] = GetUrl(chandle, str);
auto [out, res] = GetUrl(myhandle, str);
if(res != CURLE_OK) return Error(pref, MString("can't download chunk: ") + chandle.Err());
long respcode;
//curl_easy_getinfo(chandle, CURLINFO_RESPONSE_CODE, &respcode);
curl_easy_getinfo(myhandle, CURLINFO_RESPONSE_CODE, &respcode);
if(respcode == 403) out = ""; // Failed chunk download mean that this chunk contains only fill
cache->Put(str, out, 3600);
content = std::move(out);
}
if(content.Exist())
{
size_t nb, cb, bs;
blosc_cbuffer_sizes(content.Buf(), &nb, &cb, &bs);
if(cb != content.Len()) return Error(pref, MString("bytes download: ") + content.Len() + ", but compressed bytes " + cb);
if(nb != chunksize * elsize) return Error(pref, MString("decompressed bytes: ") + nb + ", but buffer size " + chunksize * elsize);
auto res = blosc_decompress_ctx(content.Buf(), data, chunksize * elsize, 1);
if(int_cast<size_t>(res) != chunksize * elsize) return Error(pref, MString("decompress only ") + res + " bytes of " + chunksize * elsize);
}
else
{
if(fill == nullptr) return Error(pref, MString("can't download chunk: ") + chandle.Err());
for(size_t i = 0; i < chunksize; i++) memcpy(michlib::P1(data) + i * elsize, fill, elsize);
}
return Error();
}
Error ZarrFunctions::Open(const MString& product, const MString& dataset, bool time)
{
static const MString pref = "Zarr::Open";
gats.clear();
dims.clear();
vars.clear();
CopernicusCatalog cat;
Json::Value json;
MString realdataset;
if(!dataset.Exist())
{
auto dsets = cat.DatasetList(product);
if(!dsets) return dsets.Add(pref, "Can't get default dataset of product " + product);
realdataset = dsets.Value()[0];
}
else
realdataset = dataset;
{
auto urlret = time ? cat.DatasetTimeURL(product, realdataset) : cat.DatasetGeoURL(product, realdataset);
if(!urlret) return urlret.Add(pref, "Can't get url for the dataset " + realdataset + " of product " + product);
url = urlret.Value();
auto ret = cat.GetJSON(url + "/.zmetadata");
if(ret)
json = ret.Value();
else
return ret.Add(pref, "can't download .zmetadata");
}
const auto& meta = json["metadata"];
if(!meta) return {pref, "No \"metadata\" key in JSON data"};
if(meta[".zattrs"]) gats = ReadAtts(meta[".zattrs"]);
auto vnames = ReadVarNames(meta);
for(size_t i = 0; i < vnames.size(); i++)
{
auto err = AddVar(vnames[i], meta[(vnames[i] + "/.zattrs").Buf()], meta[(vnames[i] + "/.zarray").Buf()]);
if(!err) return err.Add(pref, "Can't init variable " + vnames[i]);
}
return Error();
}
Loading…
Cancel
Save