diff --git a/include/MODISBINLOCAL.h b/include/MODISBINLOCAL.h new file mode 100644 index 0000000..a84bdb2 --- /dev/null +++ b/include/MODISBINLOCAL.h @@ -0,0 +1,60 @@ +#pragma once +#include "ParseArgs.h" +#include "basedata.h" +#include "mdatetime.h" +#include "mregex.h" +#include "DataAdapters/ncfilealt.h" +#include +#include +#include + +using michlib::GPL; +using michlib::MDateTime; +using michlib::RegExp; +using michlib::uint2; +using michlib::uint4; +using michlib::Floor; +using michlib::Ceil; +using michlib::pointer_cast; + +class MODISBINLOCALData +{ + MString datapath; + std::vector dataset; + std::vector times; + + bool SufFromDataset(const MString& datasetstr); + + struct Parameters: public BaseParameters + { + real lonb, latb, lone, late; + virtual ~Parameters() override = default; + }; + + public: + using Data = UngriddedData; + + MString Info() const; + MString Open(const CLArgs& args); + Data Read(const MString& var, const BaseParameters* ip, size_t tind) const; + Data ReadFile(const MString& suf, const struct Parameters* p, size_t tind) const; + + std::pair Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const; + + bool CheckVar(const MString& vname) const + { + if(dataset.empty()) return false; + bool ischl = dataset.front() == "A_CHL" || dataset.front() == "T_CHL"; + return (vname == "chl" && ischl) || (vname == "temp" && !ischl); + } + + bool isOk() const { return times.size() > 0; } + + size_t NTimes() const { return times.size(); } + + MDateTime Time(size_t i) const + { + if(!isOk() || i >= times.size()) return MDateTime(); + return times[i]; + } +}; diff --git a/include/data.h b/include/data.h index b30ea8b..0888424 100644 --- a/include/data.h +++ b/include/data.h @@ -3,10 +3,11 @@ #include "AVISOLOCAL.h" #include "BINFILE.h" #include "HYCOM.h" +#include "MODISBINLOCAL.h" #include "NEMO.h" #include -using DataVariants = std::variant; +using DataVariants = std::variant; class Data: public DataVariants { public: diff --git a/include/odm.h b/include/odm.h index 0d6c801..0802c19 100644 --- a/include/odm.h +++ b/include/odm.h @@ -30,4 +30,6 @@ constexpr bool BINFILEsup = MustSup() && MustSup static constexpr bool DisableAction = false; // Exceptions -template<> constexpr bool DisableAction = true; +template<> constexpr bool DisableAction = true; +template<> constexpr bool DisableAction = true; +template<> constexpr bool DisableAction = true; diff --git a/src/MODISBINLOCAL.cpp b/src/MODISBINLOCAL.cpp new file mode 100644 index 0000000..bad3d64 --- /dev/null +++ b/src/MODISBINLOCAL.cpp @@ -0,0 +1,262 @@ +#define MICHLIB_NOSOURCE +#include "MODISBINLOCAL.h" + +MString MODISBINLOCALData::Info() const +{ + MString ds = ""; + for(size_t i = 0; i < dataset.size(); i++) + { + if(i != 0) ds += " "; + ds += dataset[i]; + } + + // clang-format off + return + "Dataset: MODIS Level-3 Binned Data 4.6 km - " + ds + "\n" + + " Begin date: " + Time(0).ToString() + "\n" + + " End date: " + Time(NTimes()-1).ToString() + "\n" + + " Time moments: " + NTimes() + "\n" + + " Supported variables: temp, chl"; + // clang-format on +} + +MString MODISBINLOCALData::Open(const CLArgs& args) +{ + GPL.UsePrefix("MODISBINLOCAL"); + datapath = GPL.ParameterSValue("Datapath", ""); + MString var, datasetstr; + if(!args.contains("var") && !args.contains("dataset")) var = "temp"; + if(args.contains("var")) + { + var = args.at("var"); + if(var != "temp" && var != "chl") return "Unknown variable: " + var; + } + if(var.Exist() && !args.contains("dataset")) datasetstr = (var == "temp") ? "ALL" : "CHL"; + if(args.contains("dataset")) datasetstr = args.at("dataset"); + + bool ischl = SufFromDataset(datasetstr); + if(dataset.size() == 0) return "Incorrect dataset: " + datasetstr; + if(!var.Exist()) var = ischl ? "chl" : "temp"; + if((var == "temp" && ischl) || (var == "chl" && !ischl)) return "Arguments dataset and var are in contradiction"; + + MString regstring = "([0-9]{4}-[0-9]{2}-[0-9]{2})_("; + for(size_t i = 0; i < dataset.size(); i++) + { + if(i != 0) regstring += "|"; + regstring += dataset[i]; + } + regstring += ").nc"; + + RegExp regex(regstring.Buf()); + regex.Compile(); + + DIR* dir = opendir(datapath.Buf()); + struct dirent* de; + if(nullptr == dir) return "Can't open directory " + datapath; + while((de = readdir(dir))) + { + if(!regex.Match(de->d_name)) continue; + times.emplace_back(MString(de->d_name + regex.Off(1), regex.Len(1))); + } + closedir(dir); + std::sort(times.begin(), times.end()); + auto last = std::unique(times.begin(), times.end()); + times.erase(last, times.end()); + times.shrink_to_fit(); + return ""; +} + +std::pair MODISBINLOCALData::Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const +{ + std::unique_ptr ppar{new struct Parameters}; + if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"}; + ppar->lonb = args.at("lonb").ToReal(); + ppar->lone = args.at("lone").ToReal(); + ppar->latb = args.at("latb").ToReal(); + ppar->late = args.at("late").ToReal(); + + pars.SetParameter("lonb", ppar->lonb); + pars.SetParameter("latb", ppar->latb); + pars.SetParameter("lone", ppar->lone); + pars.SetParameter("late", ppar->late); + + return {ppar.release(), ""}; +} + +MODISBINLOCALData::Data MODISBINLOCALData::Read(const MString& var, const BaseParameters* ip, size_t tind) const +{ + auto p = pointer_cast(ip); + + Averager out; + + for(const auto& suf: dataset) + { + Data dat = ReadFile(suf, p, tind); + if(dat) out.Add(dat); + } + return out.Div(); +} + +MODISBINLOCALData::Data MODISBINLOCALData::ReadFile(const MString& suf, const struct MODISBINLOCALData::Parameters* p, size_t tind) const +{ + michlib::NCFileA nc; + michlib::message("Reading " + datapath + "/" + times[tind].ToString() + "_" + suf + ".nc"); + nc.Reset(datapath + "/" + times[tind].ToString() + "_" + suf + ".nc"); + if(!nc) return Data(); + MString dname; + if(suf == "A_CHL" || suf == "T_CHL") dname = "chlor_a"; + if(suf == "A_SST" || suf == "T_SST") dname = "sst"; + if(suf == "A_NSST" || suf == "T_NSST") dname = "sst"; + if(suf == "A_SST4" || suf == "T_SST4") dname = "sst4"; + + struct BinIndexType + { + uint4 start_num; + uint4 begin; + uint4 extent; + uint4 max; + }; + struct BinListType + { + uint4 bin_num; + uint2 nobs; + uint2 nscenes; + float weights; + float time_rec; + }; + struct BinDataType + { + float sum; + float sum_squared; + }; + + std::vector beg; + std::vector loc; + std::vector val; + + std::vector bins; + + // Read index + { + auto ind = nc.V("/level-3_binned_data/BinIndex"); + if(!ind) return Data(); + beg.resize(ind.DimLen(0) + 1); + beg[0] = 1; + for(size_t i = 1; i <= ind.DimLen(0); i++) beg[i] = beg[i - 1] + ind(i - 1).max; + } + + size_t ilatb, ilate; + { + size_t nr = beg.size() - 2; + ilatb = static_cast(Ceil((p->latb + 90.0) * nr / 180.0 - 0.5)); + ilate = static_cast(Floor((p->late + 90.0) * nr / 180.0 - 0.5)); + } + + for(size_t i = ilatb; i <= ilate; i++) + { + real lb = michlib::ToGeoDomainNeg(p->lonb); + real le = michlib::ToGeoDomainNeg(p->lone); + size_t nl = beg[i + 1] - beg[i]; + size_t ilonb, ilone, il; + ilonb = static_cast(Ceil((lb + 180.0) * nl / 360.0 - 0.5)); + ilone = static_cast(Floor((le + 180.0) * nl / 360.0 - 0.5)); + il = ilonb; + do bins.push_back(beg[i] + il); + while(++il % nl != ilone + 1); + } + std::sort(bins.begin(), bins.end()); + + auto InReg = [&bins = std::as_const(bins)](uint4 n) -> size_t + { + size_t a = 0, b = bins.size() - 1; + if(bins[a] == n) return a; + if(bins[b] == n) return b; + while(true) + { + size_t c = (a + b) / 2; + if(bins[c] == n) return c; + if(n < bins[c]) + b = c; + else + a = c; + if(b - a == 1) return bins.size(); + } + }; + + auto Row = [&beg = std::as_const(beg)](uint4 n) -> size_t + { + if(n >= beg.back() || n == 0) return beg.size() - 1; + size_t a = 0, b = beg.size() - 2; + if(n >= beg[a] && n < beg[a + 1]) return a; + if(n >= beg[b] && n < beg[b + 1]) return b; + while(true) + { + size_t c = (a + b) / 2; + if(n >= beg[c] && n < beg[c + 1]) return c; + if(n < beg[c]) + b = c; + else + a = c; + } + }; + auto Lat = [&Row = Row](uint4 n) + { + auto r = Row(n); + auto nr = Row(0) - 1; // Number of rows + return -90.0 + (r + 0.5) * (180.0 / nr); + }; + auto Lon = [&Row = Row, &beg = std::as_const(beg)](uint4 n) -> real + { + auto r = Row(n); + if(r == beg.size() - 1) return NAN; + return -180.0 + (n - beg[r] + 0.5) * (360.0 / (beg[r + 1] - beg[r])); + }; + + { + auto lst = nc.V("/level-3_binned_data/BinList"); + if(!lst) return Data(); + auto dat = nc.V("/level-3_binned_data/" + dname); + if(!dat) return Data(); + if(dat.DimLen(0) != lst.DimLen(0)) return Data(); + val.resize(dat.DimLen(0)); + loc.resize(lst.DimLen(0)); + for(size_t i = 0; i < dat.DimLen(0); i++) + { + loc[i] = lst(i).bin_num; + val[i] = dat(i).sum / lst(i).weights; + } + } + if(loc.size() != val.size()) return Data(); + + Data out( + bins.size(), [&bins = std::as_const(bins), &Lon = Lon](size_t i) { return Lon(bins[i]); }, [&bins = std::as_const(bins), &Lat = Lat](size_t i) { return Lat(bins[i]); }); + for(size_t i = 0; i < out.N(); i++) out(i) = out.Fillval(); + + for(size_t i = 0; i < loc.size(); i++) + { + auto ind = InReg(loc[i]); + if(ind < bins.size()) out.V(ind) = val[i]; + } + return out; +} + +bool MODISBINLOCALData::SufFromDataset(const MString& datasetstr) +{ + std::set sufs; + auto words = michlib::Split_on_words(datasetstr, ", \t\x00", false); + bool ischl = words.front() == "A_CHL" || words.front() == "T_CHL" || words.front() == "CHL"; + + for(auto ci = words.cbegin(); ci != words.cend(); ++ci) + { + if(ischl && (*ci == "A_CHL" || *ci == "CHL")) sufs.insert("A_CHL"); + if(ischl && (*ci == "T_CHL" || *ci == "CHL")) sufs.insert("T_CHL"); + if(!ischl && (*ci == "A_SST" || *ci == "SST" || *ci == "DSST" || *ci == "A_DSST" || *ci == "A_ALL" || *ci == "ALL")) sufs.insert("A_SST"); + if(!ischl && (*ci == "T_SST" || *ci == "SST" || *ci == "DSST" || *ci == "T_DSST" || *ci == "T_ALL" || *ci == "ALL")) sufs.insert("T_SST"); + if(!ischl && (*ci == "A_SST4" || *ci == "SST4" || *ci == "DSST" || *ci == "A_DSST" || *ci == "A_ALL" || *ci == "ALL")) sufs.insert("A_SST4"); + if(!ischl && (*ci == "T_SST4" || *ci == "SST4" || *ci == "DSST" || *ci == "T_DSST" || *ci == "T_ALL" || *ci == "ALL")) sufs.insert("T_SST4"); + if(!ischl && (*ci == "A_NSST" || *ci == "NSST" || *ci == "A_ALL" || *ci == "ALL")) sufs.insert("A_NSST"); + if(!ischl && (*ci == "T_NSST" || *ci == "NSST" || *ci == "T_ALL" || *ci == "ALL")) sufs.insert("T_NSST"); + } + for(const auto& w: sufs) dataset.push_back(w); + return ischl; +} diff --git a/src/data.cpp b/src/data.cpp index ba80817..679b58e 100644 --- a/src/data.cpp +++ b/src/data.cpp @@ -40,6 +40,13 @@ MString Data::Init(const CLArgs& args) if(res.Exist()) return "Can't open source " + src + ":\n" + res; *this = Data(std::move(data)); } + else if(src == "MODISBINLOCAL") + { + MODISBINLOCALData data; + auto res = data.Open(args); + if(res.Exist()) return "Can't open source " + src + ":\n" + res; + *this = Data(std::move(data)); + } else return "Unknown source: " + src; return "";