diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..ded64b9 --- /dev/null +++ b/.clang-format @@ -0,0 +1,191 @@ +--- +Language: Cpp +AccessModifierOffset: 0 +AlignAfterOpenBracket: Align +AlignArrayOfStructures: Right +AlignConsecutiveMacros: Consecutive +AlignConsecutiveAssignments: Consecutive +AlignConsecutiveBitFields: Consecutive +AlignConsecutiveDeclarations: Consecutive +AlignEscapedNewlines: Right +AlignOperands: Align +AlignTrailingComments: true +AllowAllArgumentsOnNextLine: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortEnumsOnASingleLine: true +AllowShortBlocksOnASingleLine: Always +AllowShortCaseLabelsOnASingleLine: true +AllowShortFunctionsOnASingleLine: All +AllowShortLambdasOnASingleLine: All +AllowShortIfStatementsOnASingleLine: WithoutElse +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: MultiLine +AttributeMacros: + - __capability +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterCaseLabel: true + AfterClass: true + AfterControlStatement: Always + AfterEnum: true + AfterFunction: true + AfterNamespace: true + AfterObjCDeclaration: true + AfterStruct: true + AfterUnion: true + AfterExternBlock: true + BeforeCatch: true + BeforeElse: true + BeforeLambdaBody: true + BeforeWhile: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeConceptDeclarations: true +BreakBeforeBraces: Custom +BreakBeforeInheritanceComma: false +BreakInheritanceList: AfterColon +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: AfterColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 180 +CommentPragmas: '^ IWYU pragma:' +QualifierAlignment: Leave +CompactNamespaces: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DeriveLineEnding: false +DerivePointerAlignment: false +DisableFormat: false +EmptyLineAfterAccessModifier: Never +EmptyLineBeforeAccessModifier: LogicalBlock +ExperimentalAutoDetectBinPacking: false +PackConstructorInitializers: NextLine +BasedOnStyle: '' +ConstructorInitializerAllOnOneLineOrOnePerLine: false +AllowAllConstructorInitializersOnNextLine: true +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IfMacros: + - KJ_IF_MAYBE +IncludeBlocks: Preserve +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 2 + SortPriority: 0 + CaseSensitive: false + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 3 + SortPriority: 0 + CaseSensitive: false + - Regex: '.*' + Priority: 1 + SortPriority: 0 + CaseSensitive: false +IncludeIsMainRegex: '(Test)?$' +IncludeIsMainSourceRegex: '' +IndentAccessModifiers: false +IndentCaseLabels: false +IndentCaseBlocks: false +IndentGotoLabels: false +IndentPPDirectives: None +IndentExternBlock: Indent +IndentRequires: false +IndentWidth: 1 +IndentWrappedFunctionNames: false +InsertTrailingCommas: None +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +LambdaBodyIndentation: Signature +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBinPackProtocolList: Auto +ObjCBlockIndentWidth: 2 +ObjCBreakBeforeNestedBlockParam: true +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakOpenParenthesis: 0 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 120 +PenaltyIndentedWhitespace: 0 +PointerAlignment: Left +PPIndentWidth: -1 +ReferenceAlignment: Pointer +ReflowComments: false +RemoveBracesLLVM: false +SeparateDefinitionBlocks: Leave +ShortNamespaceLines: 1 +SortIncludes: CaseSensitive +SortJavaStaticImport: Before +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: false +SpaceBeforeAssignmentOperators: true +SpaceBeforeCaseColon: false +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: false +SpaceBeforeInheritanceColon: false +SpaceBeforeParens: Never +SpaceBeforeParensOptions: + AfterControlStatements: false + AfterForeachMacros: false + AfterFunctionDefinitionName: false + AfterFunctionDeclarationName: false + AfterIfMacros: false + AfterOverloadedOperator: false + BeforeNonEmptyParentheses: false +SpaceAroundPointerQualifiers: Default +SpaceBeforeRangeBasedForLoopColon: false +SpaceInEmptyBlock: false +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: Never +SpacesInConditionalStatement: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInLineCommentPrefix: + Minimum: 1 + Maximum: -1 +SpacesInParentheses: false +SpacesInSquareBrackets: false +SpaceBeforeSquareBrackets: false +BitFieldColonSpacing: None +Standard: Latest +StatementAttributeLikeMacros: + - Q_EMIT +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TabWidth: 8 +UseCRLF: false +UseTab: Never +WhitespaceSensitiveMacros: + - STRINGIZE + - PP_STRINGIZE + - BOOST_PP_STRINGIZE + - NS_SWIFT_NAME + - CF_SWIFT_NAME +... + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..6a80759 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,92 @@ +cmake_minimum_required(VERSION 3.23) + +# Make sure the user doesn't play dirty with symlinks +get_filename_component(srcdir "${CMAKE_SOURCE_DIR}" REALPATH) +get_filename_component(bindir "${CMAKE_BINARY_DIR}" REALPATH) + +# Disallow in-source builds +if(${srcdir} STREQUAL ${bindir}) + message(FATAL_ERROR "In-source builds are not allowed." + " Please create a directory and run cmake from there, passing the path" + " to this source directory as the last argument. This process created" + " the file `CMakeCache.txt' and the directory `CMakeFiles' in ${srcdir}." + " Please remove them.") +endif() + +project(odm CXX) + +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel" FORCE) +endif() + +set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake/modules/" CACHE INTERNAL "Location of our custom CMake modules." FORCE) + +include_directories(include) + +# Checking compiler options +set(default_options -Wall) +set(linker_options -Wall) +include(CheckCXXCompilerFlag) +set(saved_CMAKE_REQUIRED_FLAGS ${CMAKE_REQUIRED_FLAGS}) +# C++20 is required +CHECK_CXX_COMPILER_FLAG(-std=c++20 COMPILER_SUPPORTS_CXX20) +if(COMPILER_SUPPORTS_CXX20) + set(default_options ${default_options} -std=c++20) +else() + message(FATAL_ERROR "The compiler ${CMAKE_CXX_COMPILER} has no C++20 support. Please use a different C++ compiler.") +endif() +# Link time optimization check +set(CMAKE_REQUIRED_FLAGS -flto) +CHECK_CXX_COMPILER_FLAG(-flto COMPILER_SUPPORTS_FLTO) +if(COMPILER_SUPPORTS_FLTO AND NOT CMAKE_BUILD_TYPE STREQUAL "Debug") + set(default_options ${default_options} -flto) + set(linker_options ${linker_options} -flto) +endif() +unset(CMAKE_REQUIRED_FLAGS) +# Clang can use lto only with gold linker +if(NOT COMPILER_SUPPORTS_FLTO) + set(CMAKE_REQUIRED_FLAGS -flto\ -fuse-ld=gold) + CHECK_CXX_COMPILER_FLAG(-flto\ -fuse-ld=gold COMPILER_SUPPORTS_FLTO_GOLD) + if(COMPILER_SUPPORTS_FLTO_GOLD AND NOT CMAKE_BUILD_TYPE STREQUAL "Debug") + set(default_options ${default_options} -flto) + set(linker_options ${linker_options} -flto\ -fuse-ld=gold) + endif() + unset(CMAKE_REQUIRED_FLAGS) +endif() +# Dwarf-4 support check +set(CMAKE_REQUIRED_FLAGS -gdwarf-4) +CHECK_CXX_COMPILER_FLAG(-gdwarf-4 COMPILER_SUPPORTS_DWARF4) +if(COMPILER_SUPPORTS_DWARF4 AND (CMAKE_BUILD_TYPE STREQUAL "Debug" OR CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo")) + set(default_options ${default_options} -gdwarf-4) + set(linker_options ${linker_options} -gdwarf-4) +endif() +unset(CMAKE_REQUIRED_FLAGS) +# -ipo flag check (Intel compiler link time optimization) +set(CMAKE_REQUIRED_FLAGS -ipo) +CHECK_CXX_COMPILER_FLAG(-ipo COMPILER_SUPPORTS_IPO) +if(COMPILER_SUPPORTS_IPO AND NOT CMAKE_BUILD_TYPE STREQUAL "Debug") + set(default_options ${default_options} -ipo) + set(linker_options ${linker_options} -ipo) +endif() +unset(CMAKE_REQUIRED_FLAGS) +# Default hidden visibility check +set(CMAKE_REQUIRED_FLAGS -fvisibility=hidden) +CHECK_CXX_COMPILER_FLAG(-fvisibility=hidden COMPILER_SUPPORTS_HIDDEN) +if(COMPILER_SUPPORTS_HIDDEN) + set(default_options ${default_options} -fvisibility=hidden) + set(linker_options ${linker_options} -fvisibility=hidden) +endif() +unset(CMAKE_REQUIRED_FLAGS) +# Export dynamic. Required as modules mast call functions from main programm +set(CMAKE_REQUIRED_FLAGS -Wl,--export-dynamic) +CHECK_CXX_COMPILER_FLAG(-Wl,--export-dynamic COMPILER_SUPPORTS_EXPORTDYNAMIC) +if(COMPILER_SUPPORTS_EXPORTDYNAMIC) + set(linker_options ${linker_options} -Wl,--export-dynamic) +else() + message(FATAL_ERROR "The compiler ${CMAKE_CXX_COMPILER} has no support of --export-dynamic. Please use a different C++ compiler.") +endif() +unset(CMAKE_REQUIRED_FLAGS) +add_compile_options(${default_options}) +set(CMAKE_REQUIRED_FLAGS ${saved_CMAKE_REQUIRED_FLAGS}) + +add_subdirectory(src) diff --git a/include/NEMO.h b/include/NEMO.h new file mode 100644 index 0000000..bcacbfa --- /dev/null +++ b/include/NEMO.h @@ -0,0 +1,312 @@ +#include "DataAdapters/ncfilealt.h" +#include "mdatetime.h" +#include "simple2ddata.h" +#include "specfunc.h" +#include "vartype.h" +#include +#include +#include + +#if !defined(M__NEMO) +#define M__NEMO + +using michlib::Ceil; +using michlib::Floor; +using michlib::int2; +using michlib::MDateTime; +using michlib::NCFileA; +using michlib::ToGeoDomainNeg; + +class NEMOData +{ + NCFileA nc, nct, ncs; + size_t xb, yb, xe, ye, layer; + std::vector depths; + std::vector times; + size_t nx; + 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; + + public: + using Data = Simple2DData; + + Data ReadVar(const NCFileA& f, const MString& name, size_t i) const + { + using DataType = int2; + constexpr DataType fill = -32767; + real unitmul = 1.0; + + auto offset = f.A(name, "add_offset"); + auto scale = f.A(name, "scale_factor"); + if(!offset || !scale) return Data(); + + auto unit = f.A(name, "units"); + if(unit && unit.Get() == "m s-1") unitmul = 100.0; + + Data data((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1), ye - yb + 1, Lonb(), Latb()); + + if(xb < xe) + { + auto var = f.V(name, {"longitude", xb, xe - xb + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); + if(!var) return Data(); + if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data(); + + for(size_t ix = 0; ix < var.DimLen(0); ix++) + for(size_t iy = 0; iy < var.DimLen(1); iy++) + { + DataType v = var(ix, iy); + data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); + } + } + else + { + auto var1 = f.V(name, {"longitude", xb}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); + auto var2 = f.V(name, {"longitude", 0, xe + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1}); + if(!(var1 && var2)) return Data(); + if((var1.DimLen(0) + var2.DimLen(0)) != data.Nx() || var1.DimLen(1) != data.Ny() || var2.DimLen(1) != data.Ny()) return Data(); + for(size_t ix = 0; ix < var1.DimLen(0); ix++) + for(size_t iy = 0; iy < var1.DimLen(1); iy++) + { + DataType v = var1(ix, iy); + data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); + } + for(size_t ix = 0; ix < var2.DimLen(0); ix++) + for(size_t iy = 0; iy < var2.DimLen(1); iy++) + { + DataType v = var2(ix, iy); + data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul); + } + } + return data; + } + + Data ReadVar(const NCFileA& f, const MString& name, const std::vector& tindex) const + { + Data out; + if(tindex.size() == 0) return out; + + std::vector count; + for(size_t i = 0; i < tindex.size(); i++) + { + Data dat = ReadVar(f, name, tindex[i]); + if(!dat) return Data(); + out.Add(dat, count); + } + out.Div(count); + return out; + } + + public: + NEMOData() = default; + + // TODO: RetVal + bool Open(const MString& type, const MString& cred, const MString& proxyurl = "") + { + if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl); + MString url, urlt, urls; + if(type == "DT") + { + url = "https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m"; + title = "NEMO Delayed time, daily mean"; + } + if(type == "NRT") + { + url = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024"; + title = "NEMO Near-real time, daily mean"; + } + if(type == "NRT6") + { + url = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-uovo"; + urlt = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-thetao"; + urls = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-so"; + title = "NEMO Near-real time, 6h resolution"; + } + nc.Reset(url); + if(!nc) return false; + + { + auto nx_ = nc.D("longitude"); + if(!nx_) + { + nc.Reset(); + return false; + } + nx = nx_; + } + + if(urlt.Exist()) + { + nct.Reset(urlt); + if(!nct) + { + nc.Reset(); + return false; + } + } + if(urls.Exist()) + { + ncs.Reset(urls); + if(!ncs) + { + nc.Reset(); + nct.Reset(); + return false; + } + } + + auto rdepths = nc.V("depth"); + if(!rdepths) + { + nc.Reset(); + nct.Reset(); + ncs.Reset(); + return false; + } + depths.resize(rdepths.DimLen(0)); + for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i); + + auto timeD = nc.V("time"); + auto timeF = nc.V("time"); + if(!(timeD || timeF)) + { + nc.Reset(); + nct.Reset(); + ncs.Reset(); + return false; + } + MDateTime refdate("1950-01-01"); + timeD ? times.resize(timeD.DimLen(0)) : times.resize(timeF.DimLen(0)); + for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast(timeD ? timeD(i) : timeF(i)) * 3600; + + return true; + } + + MString Info() const + { + MString d; + for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")"; + + // clang-format off + return + "NEMO dataset title: " + Title() + "\n" + + " Begin date: " + Time(0).ToString() + "\n" + + " End date: " + Time(NTimes()-1).ToString() + "\n" + + " Time step: " + Timestep() + " seconds\n" + + " Time moments: " + NTimes() + "\n" + + " Depths:" + d; + // clang-format on + } + + bool SetRegion(real lonbin, real lonein, real latb, real late, real depth) + { + if(!nc) return false; + real lonb = ToGeoDomainNeg(lonbin), lone = ToGeoDomainNeg(lonein); + + size_t xb_, xe_, yb_, ye_, layer_ = 0; + + yb_ = static_cast(Floor((latb + 80.0) * 12.0)); + ye_ = static_cast(Ceil((late + 80.0) * 12.0)); + if(ye_ > 2040) ye_ = 2040; + if(yb_ >= ye_) return false; + + xb_ = static_cast(Floor((lonb + 180.0) * 12.0)); + xe_ = static_cast(Ceil((lone + 180.0) * 12.0)); + if(xb_ == xe_) return false; + + if(depth < 0.0 || depth > depths.back()) + 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]) + { + layer_ = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1); + break; + } + } + xb = xb_; + xe = xe_; + yb = yb_; + ye = ye_; + layer = layer_; + return true; + } + + real Lonb() const { return nc ? (-180.0 + xb / 12.0) : -1000.0; } + real Lone() const { return nc ? (-180.0 + xe / 12.0) : -1000.0; } + real Latb() const { return nc ? (-80.0 + yb / 12.0) : -1000.0; } + real Late() const { return nc ? (-80.0 + ye / 12.0) : -1000.0; } + size_t Xb() const { return nc ? xb : 0; } + size_t Xe() const { return nc ? xe : 0; } + size_t Yb() const { return nc ? yb : 0; } + size_t Ye() const { return nc ? ye : 0; } + size_t Nx() const { return nc ? ((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1)) : 0; } + size_t Ny() const { return nc ? (ye - yb + 1) : 0; } + real Depth(size_t l) const { return nc ? depths[l] : -1000.0; } + real Depth() const { return Depth(layer); } + size_t Layer() const { return layer; } + time_t Timestep() const { return nc ? (times[1] - times[0]) : 0; } + MDateTime Time(size_t i) const + { + if((!nc) || i >= times.size()) return MDateTime(); + return times[i]; + } + + size_t NDepths() const { return depths.size(); } + size_t NTimes() const { return times.size(); } + + template Data Read(size_t it) const = delete; + + template Data Read(const std::vector& tindex) const = delete; + + MString Title() const { return title; } + + static real Fillval() { return Data::Fillval(); } +}; + +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar(nc, "uo", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar(nc, "vo", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar(nct ? nct : nc, "thetao", it); } +template<> inline NEMOData::Data NEMOData::Read(size_t it) const { return ReadVar(ncs ? ncs : nc, "so", it); } + +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar(nc, "uo", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar(nc, "vo", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar(nct ? nct : nc, "thetao", tindex); } +template<> inline NEMOData::Data NEMOData::Read(const std::vector& tindex) const { return ReadVar(ncs ? ncs : nc, "so", tindex); } + +#endif diff --git a/include/basedata.h b/include/basedata.h new file mode 100644 index 0000000..c62ee48 --- /dev/null +++ b/include/basedata.h @@ -0,0 +1,32 @@ +#include "comdefs.h" + +#if !defined(M__BASEDATA) +#define M__BASEDATA + +using michlib::real; + +class BaseData +{ + protected: + static constexpr real fillval = 1.0e10; + std::vector data; + + BaseData() = default; + + BaseData(size_t n): data(n) {} + + public: + const real& V(size_t i) const { return data[i]; } + real& V(size_t i) { return data[i]; } + + const real& operator()(size_t i) const { return data[i]; } + real& operator()(size_t i) { return data[i]; } + + size_t N() const { return data.size(); } + + static real Fillval() { return fillval; } + + explicit operator bool() const { return data.size() != 0; } +}; + +#endif diff --git a/include/odm.h b/include/odm.h new file mode 100644 index 0000000..a8c1ef2 --- /dev/null +++ b/include/odm.h @@ -0,0 +1,70 @@ +#include "BFileW.h" +#include "NEMO.h" +#include "tindexes.h" + +#if !defined(M__ODM) +#define M__ODM + +using michlib::BFileW; +using michlib::errmessage; +using michlib::GPL; +using michlib::message; +using michlib::SList; + +using DataVariants = std::variant; +using CLArgs = std::map; + +class Data: public DataVariants +{ + template static bool WriteData(BFileW& fw, const D& data) + { + for(size_t i = 0; i < data.N(); i++) + { + fw.Write(data.Lon(i)); + fw.Write(data.Lat(i)); + fw.Write(data(i) == data.Fillval() ? NAN : data(i)); + } + return true; + } + + public: + Data() = default; + + Data(DataVariants&& d): DataVariants(std::move(d)) {} + + auto Time(size_t it) const + { + return std::visit( + [it = it](const auto& arg) -> auto{ return arg.Time(it); }, *this); + } + auto NTimes() const + { + return std::visit( + [](const auto& arg) -> auto{ return arg.NTimes(); }, *this); + } + + bool Write(BFileW& fw, VarType vt, const std::vector& tindexes) const + { + return std::visit( + [&fw = fw, vt = vt, &ti = tindexes](const auto& d) -> bool + { + using T = std::decay_t; + + return std::visit( + [&fw = fw, &ti = ti, &d = d](auto v) -> bool + { + constexpr auto cvt = decltype(v)::vt; + if constexpr(isDataSupported) + return WriteData(fw, d.template Read(ti)); + else + return false; + }, + vt); + }, + *this); + } +}; + +CLArgs ParseArgs(int argc, char** argv); + +#endif diff --git a/include/simple2ddata.h b/include/simple2ddata.h new file mode 100644 index 0000000..765e775 --- /dev/null +++ b/include/simple2ddata.h @@ -0,0 +1,79 @@ +#include "basedata.h" + +#if !defined(M__SIMPLE2DDATA) +#define M__SIMPLE2DDATA + +class Simple2DData: public BaseData +{ + real x0 = 0.0, y0 = 0.0; + size_t nx = 0, ny = 0; + + public: + Simple2DData() = default; + + Simple2DData(size_t nx_, size_t ny_, real x0_, real y0_): BaseData(nx_ * ny_), x0(x0_), y0(y0_), nx(nx_), ny(ny_) {} + + const real& V(size_t i) const { return BaseData::V(i); } + real& V(size_t i) { return BaseData::V(i); } + + const real& operator()(size_t i) const { return BaseData::V(i); } + real& operator()(size_t i) { return BaseData::V(i); } + + const real& V(size_t ix, size_t iy) const { return V(iy * nx + ix); } + real& V(size_t ix, size_t iy) { return V(iy * nx + ix); } + + const real& operator()(size_t ix, size_t iy) const { return V(iy * nx + ix); } + real& operator()(size_t ix, size_t iy) { return V(iy * nx + ix); } + + size_t Nx() const { return nx; } + size_t Ny() const { return ny; } + + real Lon(size_t ix, size_t iy) const { return x0 + ix / 12.0; } + real Lat(size_t ix, size_t iy) const { return y0 + iy / 12.0; } + + real Lon(size_t i) const { return Lon(i % nx, i / nx); } + real Lat(size_t i) const { return Lat(i % nx, i / nx); } + + void InitCount(std::vector& count) const + { + count.resize(data.size()); + for(size_t i = 0; i < N(); i++) count[i] = (V(i) == fillval) ? 0 : 1; + } + + Simple2DData& Add(const Simple2DData& d, std::vector& count) + { + if(!d) return *this; + if(!*this) + { + *this = d; + InitCount(count); + } + + if(N() != d.N() || N() != count.size()) return *this; + for(size_t i = 0; i < N(); i++) + if(d.V(i) != fillval) + { + if(V(i) == fillval) + { + V(i) = d(i); + count[i] = 1; + } + else + { + V(i) += d(i); + count[i]++; + } + } + return *this; + } + + Simple2DData& Div(const std::vector& count) + { + if(N() != count.size()) return *this; + for(size_t i = 0; i < N(); i++) + if(count[i] != 0) V(i) /= count[i]; + return *this; + } +}; + +#endif diff --git a/include/tindexes.h b/include/tindexes.h new file mode 100644 index 0000000..f0d8e1f --- /dev/null +++ b/include/tindexes.h @@ -0,0 +1,67 @@ +#include "mdatetime.h" +#include +#include + +#if !defined(M__TINDEXES) +#define M__TINDEXES + +using michlib::MDateTime; + +template std::vector GetTIndexes(const T& adapter, const MDateTime& b, const MDateTime& e) +{ + std::vector out; + auto nt = adapter.NTimes(); + const MDateTime& beg = (b < e) ? b : e; + const MDateTime& end = (b > e) ? b : e; + + if(beg > adapter.Time(nt - 1) || end < adapter.Time(0)) return out; + + size_t ib = 0, ie = nt - 1; + for(size_t i = 0; i < nt; i++) + if(adapter.Time(i) >= beg) + { + ib = i; + break; + } + + for(size_t i = nt; i != 0; i--) + if(adapter.Time(i - 1) <= end) + { + ie = i - 1; + break; + } + + out.resize(ie - ib + 1); + for(size_t i = 0; i < ie - ib + 1; i++) out[i] = i + ib; + return out; +} + +template std::vector GetTIndexes(const T& adapter, const MString& regex) +{ + std::vector out; + + std::unique_ptr regbuf(new regex_t); + if(0 != regcomp(regbuf.get(), regex.Buf(), REG_EXTENDED | REG_NOSUB)) return out; + + for(size_t i = 0; i < adapter.NTimes(); i++) + { + MString date = adapter.Time(i).ToString(); + if(0 != regexec(regbuf.get(), date.Buf(), 0, 0, 0)) continue; + out.push_back(i); + } + + return out; +} + +template size_t GetTIndex(const T& adapter, const MDateTime& t) +{ + size_t nt = adapter.NTimes(); + + if(t <= adapter.Time(0)) return 0; + if(t >= adapter.Time(nt - 1)) return nt - 1; + for(size_t i = 0; i < nt - 1; i++) + if(t >= adapter.Time(i) && t <= adapter.Time(i + 1)) return (t - adapter.Time(i) <= adapter.Time(i + 1) - t) ? i : (i + 1); + return 0; +} + +#endif diff --git a/include/vartype.h b/include/vartype.h new file mode 100644 index 0000000..6232973 --- /dev/null +++ b/include/vartype.h @@ -0,0 +1,96 @@ +#include "MString.h" +#include + +#if !defined(M__VARTYPE) +#define M__VARTYPE + +using michlib::MString; + +namespace vartype +{ +enum class Vartype +{ + NONE, + U, + V, + TEMP, + SAL, + CHL +}; + +template struct VartypeCarrier +{ + static const Vartype vt = vt_; +}; + +using VartypeUnion = std::variant, VartypeCarrier, VartypeCarrier, VartypeCarrier, + VartypeCarrier, VartypeCarrier>; + +template +concept HasVar = requires(T t) +{ + { + t.template Read(0)(0, 0) + } -> std::convertible_to; + { + t.template Read(std::vector())(0, 0) + } -> std::convertible_to; +}; +} // namespace vartype + +template static constexpr bool isDataSupported = vartype::HasVar; + +class VarType: public vartype::VartypeUnion +{ + auto VT() const + { + return std::visit( + [](auto v) -> auto{ return decltype(v)::vt; }, *this); + } + + public: + template VarType(vartype::VartypeCarrier&& vc): vartype::VartypeUnion(std::move(vc)) {} + + VarType(MString str) + { + str.ToLower(); + if(str == "u") *this = vartype::VartypeCarrier(); + if(str == "v") *this = vartype::VartypeCarrier(); + if(str == "t" || str == "temp" || str == "temperature") *this = vartype::VartypeCarrier(); + if(str == "s" || str == "sal" || str == "salinity") *this = vartype::VartypeCarrier(); + if(str == "c" || str == "chl" || str == "chlorophyll") *this = vartype::VartypeCarrier(); + } + + MString Name() const + { + switch(VT()) + { + case(vartype::Vartype::NONE): return "none"; + case(vartype::Vartype::U): return "U"; + case(vartype::Vartype::V): return "V"; + case(vartype::Vartype::TEMP): return "Temperature"; + case(vartype::Vartype::SAL): return "Salinity"; + case(vartype::Vartype::CHL): return "Chlorophyll"; + } + return "none"; + } + + bool Ok() const { return VT() != vartype::Vartype::NONE; } + explicit operator bool() const { return Ok(); } + bool operator==(const VarType& vt) const { return VT() == vt.VT(); } + bool operator!=(const VarType& vt) const { return VT() != vt.VT(); } + template bool isSupported() const + { + return std::visit( + [](const auto v) -> bool + { + if constexpr(isDataSupported) + return true; + else + return false; + }, + *this); + } +}; + +#endif diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..f8a830c --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,10 @@ +set(EXENAME odm) + +find_library(netcdf netcdf REQUIRED) + +include_directories(../michlib/michlib) + +file(GLOB srcs *.cpp) +add_executable(${EXENAME} ${srcs}) +target_link_libraries(${EXENAME} ${linker_options} ${netcdf}) +set_target_properties(${EXENAME} PROPERTIES POSITION_INDEPENDENT_CODE ON) diff --git a/src/ParseArgs.cpp b/src/ParseArgs.cpp new file mode 100644 index 0000000..206043e --- /dev/null +++ b/src/ParseArgs.cpp @@ -0,0 +1,29 @@ +#define MICHLIB_NOSOURCE +#include "GPL.h" +#include "odm.h" + +CLArgs ParseArgs(int argc, char** argv) +{ + CLArgs out; + + out["_name"] = argv[0]; + for(int i = 0; i < argc; i++) + { + MString carg = argv[i]; + out["_arg" + MString(i)] = carg; + + for(size_t c = 0; c < carg.Len(); c++) + if(carg[c] == '=') + { + MString name = carg.SubStr(1, c); + MString val = carg.SubStr(c + 2, carg.Len() - c - 1); + out[name] = val; + break; + } + } + + SList sl; + michlib_internal::ParseParameterFile("/etc/michlibconf", sl, 0); + + return out; +} diff --git a/src/odm.cpp b/src/odm.cpp new file mode 100644 index 0000000..068ef66 --- /dev/null +++ b/src/odm.cpp @@ -0,0 +1,168 @@ +#include "odm.h" +#include "GPL.h" + +inline void Usage(const MString& arg0) +{ + message(arg0 + " (=...)"); + message("Keys are:"); + message(" source. Required. May be: NEMO"); + message(" var. Required. May be: temp, sal or chl"); + message(" time. Time moment or regular expression. If present, timeb and timee must be absent"); + message(" timeb, timee. Time interval. If present, time must be absent"); + message(" out. Output file. Default: out.bin"); + message(" Keys for source=NEMO"); + message(" var can be only temp or sal"); + message(" dataset. Can be DT, NRT or NRT6. Default: DT"); + message(" cred. Login and password for COPERNICUS site. Default works"); + message(" proxy. Proxy for access to COPERNICUS. Default works"); + message(" layer and/or depth. Layer or depth of NEMO dataset. If depth is specified, layer is ignored. Default: layer=0"); + message(" lonb, lone, latb, late. Required. Region of interest"); +} + +int main(int argc, char** argv) +{ + if(argc == 1) + { + Usage(argv[0]); + return 2; + } + auto args = ParseArgs(argc, argv); + + if(!args.contains("source")) + { + errmessage("No source specified!"); + return 1; + } + + VarType vtype(args["var"]); + if(!vtype) + { + errmessage("Incorrect or no variable specified"); + return 1; + } + + michlib_internal::ParameterListEx pars; + pars.UsePrefix(""); + + pars.SetParameter("source", args["source"]); + pars.SetParameter("variable", vtype.Name()); + + Data data; + + if(args["source"] == "NEMO") + { + if(!vtype.isSupported()) + { + errmessage("Variable " + args["var"] + " is unsupported by NEMO"); + return 1; + } + + GPL.UsePrefix("AVISO"); + NEMOData ndata; + MString dataset = args["dataset"]; + if(!dataset.Exist()) dataset = "DT"; + MString cred = args["cred"]; + if(!cred.Exist()) cred = GPL.ParameterSValue("COPERNICUS_USER", ""); + MString proxy; + if(args.contains("proxy")) + MString proxy = args["proxy"]; + else + proxy = GPL.ParameterSValue("COPERNICUS_PROXY", ""); + + if(!ndata.Open(dataset, cred, proxy)) + { + errmessage("Can't open NEMO dataset"); + return 1; + } + + size_t layer = 0; + if(args.contains("layer")) layer = args["layer"].ToInteger(); + if(!args.contains("depth") && layer >= ndata.NDepths()) + { + errmessage(MString("Layer ") + layer + " is too deep!"); + return 1; + } + real depth = args.contains("depth") ? args["depth"].ToReal() : ndata.Depth(layer); + + if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) + { + errmessage("Region not specified (lonb, lone, latb, late)"); + return 1; + } + if(!ndata.SetRegion(args["lonb"].ToReal(), args["lone"].ToReal(), args["latb"].ToReal(), args["late"].ToReal(), depth)) + { + errmessage("Can't set region"); + return 1; + } + + pars.SetParameter("depth", depth); + pars.SetParameter("layer", ndata.Layer()); + pars.SetParameter("dataset", dataset); + pars.SetParameter("lonb", args["lonb"].ToReal()); + pars.SetParameter("latb", args["latb"].ToReal()); + pars.SetParameter("lone", args["lone"].ToReal()); + pars.SetParameter("late", args["late"].ToReal()); + + data = DataVariants(std::move(ndata)); + } + else + { + errmessage("Unknown source " + args["source"]); + return 1; + } + + if(args.contains("time") && (args.contains("timeb") || args.contains("timee"))) + { + errmessage("Time must be set via time parameter or timeb and timee parameter but not via both"); + return 1; + } + + if(!(args.contains("time") || (args.contains("timeb") && args.contains("timee")))) + { + errmessage("Time must be set via time parameter or timeb and timee parameter"); + return 1; + } + + BFileW fw; + MString name = args["out"]; + if(!name.Exist()) name = "out.bin"; + + fw.Create(name, 3); + /* + if(!fw.Create(name,3)) + { + errmessage("Can't create file "+name); + return 1; + } + */ + fw.SetParameters(pars); + fw.UsePrefix(""); + + if(args.contains("time")) + { + MDateTime time; + if(time.FromString(args["time"])) // One date + { + auto it = GetTIndex(data, time); + fw.SetParameter("time", data.Time(it).ToString()); + data.Write(fw, vtype, std::vector(1, it)); + } + else // Regular expression + { + fw.SetParameter("timeregex", args["time"]); + data.Write(fw, vtype, GetTIndexes(data, args["time"])); + } + } + else // Bdate, edate + { + MDateTime tb(args["timeb"]), te(args["timee"]); + fw.SetParameter("timeb", tb.ToString()); + fw.SetParameter("timee", te.ToString()); + data.Write(fw, vtype, GetTIndexes(data, tb, te)); + } + + fw.Finalize(); + fw.Close(); + + return 0; +}