From 467a24e51ce4fcebb5753e10aced05689a0ecd59 Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Thu, 20 Jul 2023 15:54:52 +1000 Subject: [PATCH] Support for units of measurement in sources and the tsc module. --- actions/actiontsc.h | 5 +++-- include/basedata.h | 9 +++++++-- include/simple2ddata.h | 5 ++++- sources/AVISOLOCAL.cpp | 15 ++++++++++++--- sources/BINFILE.cpp | 3 ++- sources/MODISBINLOCAL.cpp | 7 ++++++- src/layereddata.cpp | 18 +++++++++++++++--- 7 files changed, 49 insertions(+), 13 deletions(-) diff --git a/actions/actiontsc.h b/actions/actiontsc.h index 1bf97f3..6726a97 100644 --- a/actions/actiontsc.h +++ b/actions/actiontsc.h @@ -1,6 +1,6 @@ #pragma once -#include "actiondep.h" #include "BFileW.h" +#include "actiondep.h" #include using michlib::BFileW; @@ -47,6 +47,7 @@ template MString ActionTSC::DoAction(const CLArgs& args, D& ds) auto data = Read(ds, vname, p, tindexes); if(!data) return "Can't read data"; + if(!data.Unit().Exist()) michlib::errmessage("Unknown measurement unit!"); BFileW fw; MString name = args.contains("out") ? args.at("out") : ""; @@ -55,7 +56,7 @@ template MString ActionTSC::DoAction(const CLArgs& args, D& ds) fw.Create(name, 3); fw.SetColumnName(1, "Longitude"); fw.SetColumnName(2, "Latitude"); - fw.SetColumnName(3, vname); + fw.SetColumnName(3, vname + ", " + (data.Unit().Exist() ? data.Unit() : "unknown")); fw.SetParameters(pars); for(size_t i = 0; i < data.N(); i++) { diff --git a/include/basedata.h b/include/basedata.h index f1544cf..922d3f6 100644 --- a/include/basedata.h +++ b/include/basedata.h @@ -20,8 +20,9 @@ class BaseData protected: static constexpr real fillval = 1.0e10; std::vector data; + MString unit; - BaseData(size_t n): data(n) {} + BaseData(size_t n, MString&& unit_ = ""): data(n), unit(std::move(unit_)) {} public: BaseData() = default; @@ -39,6 +40,10 @@ class BaseData static real Fillval() { return fillval; } explicit operator bool() const { return N() != 0; } + + void SetUnit(const MString& str) { unit = str; } + + const MString& Unit() const { return unit; } }; class UngriddedData: public BaseData @@ -46,7 +51,7 @@ class UngriddedData: public BaseData std::vector lons, lats; public: - template UngriddedData(size_t n, Lon genlon, Lat genlat): BaseData(n), lons(n), lats(n) + template UngriddedData(size_t n, Lon genlon, Lat genlat, MString&& unit = ""): BaseData(n, std::move(unit)), lons(n), lats(n) { for(size_t i = 0; i < n; i++) { diff --git a/include/simple2ddata.h b/include/simple2ddata.h index d24ef6e..80bc958 100644 --- a/include/simple2ddata.h +++ b/include/simple2ddata.h @@ -10,7 +10,10 @@ class Simple2DData: public BaseData public: Simple2DData() = default; - Simple2DData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_): BaseData(nx_ * ny_), x0(x0_), y0(y0_), nx(nx_), ny(ny_), xstep(xs_), ystep(ys_) {} + Simple2DData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_, MString&& unit = ""): + BaseData(nx_ * ny_, std::move(unit)), x0(x0_), y0(y0_), nx(nx_), ny(ny_), xstep(xs_), ystep(ys_) + { + } const real& V(size_t i) const { return BaseData::V(i); } real& V(size_t i) { return BaseData::V(i); } diff --git a/sources/AVISOLOCAL.cpp b/sources/AVISOLOCAL.cpp index 79e6f7c..cc27f5d 100644 --- a/sources/AVISOLOCAL.cpp +++ b/sources/AVISOLOCAL.cpp @@ -123,6 +123,7 @@ AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParame auto v = Read("v", ip, i); if(!(u && v)) return Data(); auto out = u; + out.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit()); for(size_t ind = 0; ind < out.N(); ind++) { if(u.IsFill(ind) || v.IsFill(ind)) @@ -196,10 +197,18 @@ template AVISOLOCALData::Data AVISOLOCALData::ReadVarRaw(const N if(xb == xe) return Data(); - auto unit = nc.A(name, "units"); - if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0; + MString unit; + { + auto unitatt = nc.A(name, "units"); + if(unitatt) unit = unitatt; + } + if(unit == "m s-1" || unit == "m/s") + { + unitmul = 100.0; + unit = "cm/s"; + } - Data data((xb < xe) ? (xe - xb + 1) : (dn.nx + xe - xb + 1), ye - yb + 1, lons(xb), lats(yb), lonstep, latstep); + Data data((xb < xe) ? (xe - xb + 1) : (dn.nx + xe - xb + 1), ye - yb + 1, lons(xb), lats(yb), lonstep, latstep, std::move(unit)); if(xb < xe) { diff --git a/sources/BINFILE.cpp b/sources/BINFILE.cpp index 92a0f1b..9b628a9 100644 --- a/sources/BINFILE.cpp +++ b/sources/BINFILE.cpp @@ -60,7 +60,7 @@ BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const // Only rectangular grids are supported real xs = (data->Lon(data->Nx() - 1, data->Ny() - 1) - data->Lon(0, 0)) / (data->Nx() - 1); real ys = (data->Lat(data->Nx() - 1, data->Ny() - 1) - data->Lat(0, 0)) / (data->Ny() - 1); - Data out(data->Nx(), data->Ny(), data->Lon(0, 0), data->Lat(0, 0), xs, ys); + Data out(data->Nx(), data->Ny(), data->Lon(0, 0), data->Lat(0, 0), xs, ys, "cm/s"); // TODO: Information about unit must be taken from datafile // U and U2 from u and v if(vname == "U" || vname == "U2") @@ -70,6 +70,7 @@ BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const auto v = Read("v", i); if(!(u && v)) return Data(); auto out = u; + out.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit()); for(size_t ind = 0; ind < out.N(); ind++) { if(u.IsFill(ind) || v.IsFill(ind)) diff --git a/sources/MODISBINLOCAL.cpp b/sources/MODISBINLOCAL.cpp index a6e5f60..dd09451 100644 --- a/sources/MODISBINLOCAL.cpp +++ b/sources/MODISBINLOCAL.cpp @@ -103,11 +103,15 @@ MODISBINLOCALData::Data MODISBINLOCALData::ReadFile(const MString& suf, const st michlib::message("Reading " + datapath + "/" + times[tind].ToString() + "_" + suf + ".nc"); nc.Reset(datapath + "/" + times[tind].ToString() + "_" + suf + ".nc"); if(!nc) return Data(); - MString dname; + MString dname, unit; 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"; + if(suf == "A_CHL" || suf == "T_CHL") + unit = "mg m-3"; + else + unit = "degrees_C"; struct BinIndexType { @@ -231,6 +235,7 @@ MODISBINLOCALData::Data MODISBINLOCALData::ReadFile(const MString& suf, const st 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(); + out.SetUnit(unit); for(size_t i = 0; i < loc.size(); i++) { diff --git a/src/layereddata.cpp b/src/layereddata.cpp index 1ddfc28..8ea5fa1 100644 --- a/src/layereddata.cpp +++ b/src/layereddata.cpp @@ -194,6 +194,7 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters* auto sal = Read("sal", ip, i); if(!(temp && sal)) return Data(); auto out = temp; + out.SetUnit("degrees_C"); for(size_t ind = 0; ind < out.N(); ind++) { if(temp.IsFill(ind) || sal.IsFill(ind)) @@ -210,6 +211,7 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters* auto sal = Read("sal", ip, i); if(!(temp && sal)) return Data(); auto out = temp; + out.SetUnit("degrees_C"); for(size_t ind = 0; ind < out.N(); ind++) { if(temp.IsFill(ind) || sal.IsFill(ind)) @@ -227,6 +229,7 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters* auto sal = Read("sal", ip, i); if(!(temp && sal)) return Data(); auto out = temp; + out.SetUnit("kg m-3"); for(size_t ind = 0; ind < out.N(); ind++) { if(temp.IsFill(ind) || sal.IsFill(ind)) @@ -245,6 +248,7 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters* auto v = Read("v", ip, i); if(!(u && v)) return Data(); auto out = u; + out.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit()); for(size_t ind = 0; ind < out.N(); ind++) { if(u.IsFill(ind) || v.IsFill(ind)) @@ -299,10 +303,18 @@ template LayeredData::Data LayeredData::ReadVarRaw(const NC& f, if(a_scale_f) scale = a_scale_f; } - auto unit = f->A(name, "units"); - if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0; + MString unit; + { + auto unitatt = f->A(name, "units"); + if(unitatt) unit = unitatt; + } + 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); + 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)); if(p->xb < p->xe) {