From 09f13ac243516c55e69f309b53d8d4d7699de03a Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Tue, 5 Sep 2023 16:01:31 +1000 Subject: [PATCH] Gradient calculation code (without filtering) --- actions/actiongrad.cpp | 219 +++++++++++++++++++++++++++++++++++++++++ actions/actiongrad.h | 146 +++++++++++++++++++++++++++ sources/TSCDATA.cpp | 129 ++++++++++++++++++++++++ sources/TSCDATA.h | 35 +++++++ 4 files changed, 529 insertions(+) create mode 100644 actions/actiongrad.cpp create mode 100644 actions/actiongrad.h create mode 100644 sources/TSCDATA.cpp create mode 100644 sources/TSCDATA.h diff --git a/actions/actiongrad.cpp b/actions/actiongrad.cpp new file mode 100644 index 0000000..1b69be0 --- /dev/null +++ b/actions/actiongrad.cpp @@ -0,0 +1,219 @@ +#define MICHLIB_NOSOURCE +#include "actiongrad.h" + +MString GradMethods::NCFileW::Create(const MString& name, const MString& history, const std::vector& vnames, const std::vector& lnames, + const std::vector& lons, const std::vector& lats, int compress) +{ + int ret; + + const float node_offset = 0.0; + + auto nx = lons.size(); + auto ny = lats.size(); + + int xid, yid; + + MString text; + + struct + { + union + { + struct + { + int ydimid, xdimid; + }; + int dimid[2]; + }; + } dim; + + struct Keeper + { + bool active; + int ncid; + ~Keeper() + { + if(active) nc_close(ncid); + } + operator int() { return ncid; } + } newnc; + + const auto fill = static_cast(std::numeric_limits::max()); + + // Creating + ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &newnc.ncid); + if(ret != NC_NOERR) return "Can't create netcdf file: " + name; + newnc.active = true; + + // Global attributes + ret = nc_put_att_text(newnc, NC_GLOBAL, "history", history.Len() + 1, history.Buf()); + if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file"; + + ret = nc_put_att(newnc, NC_GLOBAL, "node_offset", NC_FLOAT, 1, &node_offset); + if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file"; + + // Dimensions + ret = nc_def_dim(newnc, "longitude", nx, &dim.xdimid); + if(ret != NC_NOERR) return "Can't create x-dimension in the netcdf file"; + ret = nc_def_dim(newnc, "latitude", ny, &dim.ydimid); + if(ret != NC_NOERR) return "Can't create y-dimension in the netcdf file"; + + // lon, lat variables + ret = nc_def_var(newnc, "longitude", NC_FLOAT, 1, &dim.xdimid, &xid); + if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file"; + ret = nc_def_var(newnc, "latitude", NC_FLOAT, 1, &dim.ydimid, &yid); + if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file"; + + ret = nc_def_var_deflate(newnc, xid, 1, 1, compress); + if(ret != NC_NOERR) return "Can't set deflate parameters for longitude variable in the netcdf file"; + ret = nc_def_var_deflate(newnc, yid, 1, 1, compress); + if(ret != NC_NOERR) return "Can't set deflate parameters for latitude variable in the netcdf file"; + + text = "longitude"; + ret = nc_put_att_text(newnc, xid, "standard_name", text.Len() + 1, text.Buf()); + if(ret != NC_NOERR) return "Can't write standard_name attribute of longitude variable in the netcdf file"; + text = "latitude"; + ret = nc_put_att_text(newnc, yid, "standard_name", text.Len() + 1, text.Buf()); + if(ret != NC_NOERR) return "Can't write standard_name attribute of latitude variable in the netcdf file"; + + text = "Longitude"; + ret = nc_put_att_text(newnc, xid, "long_name", text.Len() + 1, text.Buf()); + if(ret != NC_NOERR) return "Can't write long_name attribute of longitude variable in the netcdf file"; + text = "Latitude"; + ret = nc_put_att_text(newnc, yid, "long_name", text.Len() + 1, text.Buf()); + if(ret != NC_NOERR) return "Can't write long_name attribute of latitude variable in the netcdf file"; + + // Variables + for(size_t i = 0; i < vnames.size(); i++) + { + int vid; + ret = nc_def_var(newnc, vnames[i].Buf(), NC_FLOAT, 2, dim.dimid, &vid); + if(ret != NC_NOERR) return "Can't create " + vnames[i] + " variable in the netcdf file"; + + ret = nc_def_var_deflate(newnc, vid, 1, 1, compress); + if(ret != NC_NOERR) return "Can't set deflate parameters for " + vnames[i] + " variable in the netcdf file"; + + if(lnames[i].Exist()) ret = nc_put_att_text(newnc, vid, "long_name", lnames[i].Len() + 1, lnames[i].Buf()); + if(ret != NC_NOERR) return "Can't write long_name attribute of " + vnames[i] + " variable in the netcdf file"; + + ret = nc_put_att_float(newnc, vid, "_FillValue", NC_FLOAT, 1, &fill); + if(ret != NC_NOERR) return "Can't write _FillValue attribute of " + vnames[i] + " variable in the netcdf file"; + } + + // End definitions + nc_enddef(newnc); + + // Writing lon, lat + const size_t i = 0; + + ret = nc_put_vara_float(newnc, xid, &i, &nx, lons.data()); + if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file"; + ret = nc_put_vara_float(newnc, yid, &i, &ny, lats.data()); + if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file"; + + ncid = newnc; + newnc.active = false; + return ""; +} + +MString GradMethods::NCFileW::WriteVar(const MString& name, const std::vector& data) +{ + int ret; + int vid; + int xid, yid; + size_t nx, ny; + + ret = nc_inq_varid(ncid, name.Buf(), &vid); + if(ret != NC_NOERR) return "Can't find variable " + name + " in the netcdf file"; + + ret = nc_inq_dimid(ncid, "longitude", &xid); + if(ret != NC_NOERR) return "Can't find longitude dimension in the netcdf file"; + ret = nc_inq_dimid(ncid, "latitude", &yid); + if(ret != NC_NOERR) return "Can't find latitude dimension in the netcdf file"; + + ret = nc_inq_dimlen(ncid, xid, &nx); + if(ret != NC_NOERR) return "Can't find longitude size in the netcdf file"; + ret = nc_inq_dimlen(ncid, yid, &ny); + if(ret != NC_NOERR) return "Can't find latitude size in the netcdf file"; + + const size_t i[2] = {0, 0}; + const size_t c[2] = {ny, nx}; + + GradMethods::DataType buf[nx * ny]; + for(size_t i = 0; i < nx * ny; i++) buf[i] = static_cast(data[i]); + + ret = nc_put_vara_float(ncid, vid, i, c, buf); + if(ret != NC_NOERR) return "Can't write " + name + " variable in the netcdf file"; + return ""; +} + +GradMethods::Matrix::Matrix(const std::vector& in, size_t nx_, size_t ny_, struct GradMethods::MinMax minmax): nx(nx_), ny(ny_), data(nx_ * ny_) +{ + if(minmax.automin || minmax.automax) + { + DataType min = in[0]; + DataType max = in[0]; + for(size_t i = 1; i < in.size(); i++) + if(in[i] != minmax.fill) + { + min = std::min(min, in[i]); + max = std::max(max, in[i]); + } + if(minmax.automin) minmax.min = min; + if(minmax.automax) minmax.max = max; + } + + if(minmax.log) + { + minmax.min = michlib_internal::RealType::Log(minmax.min); + minmax.max = michlib_internal::RealType::Log(minmax.max); + } + + DataType a = (std::numeric_limits::max() - 1) / (minmax.max - minmax.min); + for(size_t i = 1; i < in.size(); i++) + { + DataType v = minmax.log ? michlib_internal::RealType::Log(in[i]) : in[i]; + if(in[i] == minmax.fill) + data[i] = std::numeric_limits::max(); + else if(v <= minmax.min) + data[i] = 0; + else if(v >= minmax.max) + data[i] = std::numeric_limits::max() - 1; + else + data[i] = static_cast(michlib_internal::RealType::Round(a * (v - minmax.min))); + } +} + +void GradMethods::Matrix::Grad() +{ + std::vector out(data.size()); + const auto bad = std::numeric_limits::max(); + + for(size_t iy = 0; iy < ny; iy++) + for(size_t ix = 0; ix < nx; ix++) + { + if(iy < 1 || ix < 1 || iy > ny - 2 || ix > nx - 2) + out[iy * nx + ix] = bad; + else if(V(ix - 1, iy - 1) == bad || V(ix, iy - 1) == bad || V(ix + 1, iy - 1) == bad || V(ix - 1, iy) == bad || V(ix, iy) == bad || V(ix + 1, iy) == bad || + V(ix - 1, iy + 1) == bad || V(ix, iy + 1) == bad || V(ix + 1, iy + 1) == bad) + out[iy * nx + ix] = bad; + else + { + using IT = michlib::int4; + // Possible but unlikely overflow + const IT m1 = -1; + const IT m2 = -2; + const IT p1 = 1; + const IT p2 = 2; + + IT gx = m1 * V(ix - 1, iy + 1) + p1 * V(ix + 1, iy + 1) + m2 * V(ix - 1, iy) + p2 * V(ix + 1, iy) + m1 * V(ix - 1, iy - 1) + p1 * V(ix + 1, iy - 1); + IT gy = m1 * V(ix - 1, iy - 1) + p1 * V(ix - 1, iy + 1) + m2 * V(ix, iy - 1) + p2 * V(ix, iy + 1) + m1 * V(ix + 1, iy - 1) + p1 * V(ix + 1, iy + 1); + + auto sq = static_cast(michlib::Round(michlib::Hypot(gx, gy))); + if(sq >= bad) sq = bad - 1; + out[iy * nx + ix] = static_cast(sq); + } + } + + data = std::move(out); +} diff --git a/actions/actiongrad.h b/actions/actiongrad.h new file mode 100644 index 0000000..64edc89 --- /dev/null +++ b/actions/actiongrad.h @@ -0,0 +1,146 @@ +#pragma once +#include "actiondep.h" +#include "ncfilew.h" +#include "ncfuncs.h" + +class GradMethods +{ + public: + using DataType = float; + + static DataType ToNum(const MString& str) { return michlib_internal::RealType::String2Real(str.Buf()); } + + struct MinMax + { + bool automin, automax, log; + DataType min, max; + DataType fill; + }; + + class Matrix; + class NCFileW; +}; + +class GradMethods::Matrix +{ + public: + using MDataType = michlib::uint2; + + private: + size_t nx, ny; + std::vector data; + + public: + Matrix(const std::vector& in, size_t nx_, size_t ny_, struct MinMax minmax); + + void Grad(); + + auto Nx() const { return nx; } + auto Ny() const { return ny; } + + const auto& V(size_t ix, size_t iy) const { return data[iy * nx + ix]; } + + auto& V(size_t ix, size_t iy) { return data[iy * nx + ix]; } + + const auto& Data() const { return data; } +}; + +class GradMethods::NCFileW +{ + bool opened; + int ncid; + + public: + NCFileW(): opened(false) {} + + MString Create(const MString& name, const MString& history, const std::vector& vnames, const std::vector& lnames, const std::vector& lons, + const std::vector& lats, int compress); + MString WriteVar(const MString& name, const std::vector& data); + + ~NCFileW() { Close(); } + + void Close() + { + if(opened) nc_close(ncid); + opened = false; + } +}; + +template +concept GradSupported = requires(T t, const MString& vname) { + { + t.ReadVar(vname) + } -> std::same_as>; +}; + +ADD_ACTION(GRAD, grad, GradSupported, GradMethods); + +template MString ActionGRAD::DoAction(const CLArgs& args, D& ds) +{ + auto resop = ds.Open(args); + if(resop.Exist()) return "Can't open source: " + resop; + + MString name = args.contains("out") ? args.at("out") : "out.nc"; + + MString min = args.contains("min") ? args.at("min") : "auto"; + MString max = args.contains("max") ? args.at("max") : "auto"; + + int compress = args.contains("compress") ? args.at("compress").ToInt() : 3; + + std::vector data; + + std::vector lnames; + + // Read data + for(size_t i = 0; i < ds.NVar(); i++) + { + const MString& name = ds.VarNames()[i]; + const MString& lname = ds.LongNames()[i]; + bool hmin = args.contains(name + "_min"); + bool hmax = args.contains(name + "_max"); + struct MinMax minmax; + + minmax.log = args.contains(name + "_log"); + + if(hmin) + { + MString vmin = args.at(name + "_min"); + minmax.automin = (vmin == "auto"); + minmax.min = ToNum(vmin); + } + else + { + minmax.automin = (min == "auto"); + minmax.min = ToNum(min); + } + + if(hmax) + { + MString vmax = args.at(name + "_max"); + minmax.automax = (vmax == "auto"); + minmax.max = ToNum(vmax); + } + else + { + minmax.automax = (max == "auto"); + minmax.max = ToNum(max); + } + + minmax.fill = ds.FillVal(name); + + data.emplace_back(ds.ReadVar(name), ds.Nx(), ds.Ny(), minmax); + lnames.emplace_back(lname + ", gradient"); + } + + NCFileW fw; + fw.Create(name, (ds.History().Exist() ? (ds.History() + "; ") : "") + args.at("_cmdline"), ds.VarNames(), lnames, ds.ReadLons(), ds.ReadLats(), compress); + + for(size_t i = 0; i < ds.NVar(); i++) + { + const MString& name = ds.VarNames()[i]; + data[i].Grad(); + fw.WriteVar(name, data[i].Data()); + } + + return ""; +}; diff --git a/sources/TSCDATA.cpp b/sources/TSCDATA.cpp new file mode 100644 index 0000000..1326f7d --- /dev/null +++ b/sources/TSCDATA.cpp @@ -0,0 +1,129 @@ +#define MICHLIB_NOSOURCE +#include "TSCDATA.h" + +MString TSCDATAData::Open(const CLArgs& args) +{ + if(!args.contains("dataset")) return "path to data not specified"; + MString dataset = args.at("dataset"); + + michlib::NCFileA newnc; + std::vector newvnames, newlnames; + MString newhistory; + + newnc.Reset(dataset); + if(!newnc) return "Can't open file " + dataset; + + auto head = newnc.Header(); + if(head.Dimensions().size() != 2) return "Unsupported number of dimensions"; + if((head.Dimensions()[0].Name() != "longitude" || head.Dimensions()[1].Name() != "latitude") && + (head.Dimensions()[1].Name() != "longitude" || head.Dimensions()[0].Name() != "latitude")) + return "Unsupported dimensions names"; + + if(head.Dimensions()[0].Name() == "longitude") + { + nx = head.Dimensions()[0].Len(); + ny = head.Dimensions()[1].Len(); + } + else + { + ny = head.Dimensions()[0].Len(); + nx = head.Dimensions()[1].Len(); + } + + { + bool lonfound = false, latfound = false; + for(const auto& v: head.Variables()) + if(v.Dimensions().size() == 1 && v.Type().Id() == NC_FLOAT) + { + lonfound = lonfound || (v.Dimensions()[0].Name() == "longitude" && v.Name() == "longitude"); + latfound = latfound || (v.Dimensions()[0].Name() == "latitude" && v.Name() == "latitude"); + } + if(!lonfound) return "Longitude not found"; + if(!latfound) return "Latitude not found"; + } + + for(const auto& v: head.Variables()) + { + if(v.Dimensions().size() != 2) continue; + if(v.Type().Id() != NC_FLOAT) continue; + if((v.Dimensions()[0].Name() != "longitude" || v.Dimensions()[1].Name() != "latitude") && (v.Dimensions()[1].Name() != "longitude" || v.Dimensions()[0].Name() != "latitude")) + continue; + newvnames.push_back(v.Name()); + auto lname = newnc.A(v.Name(), "long_name"); + newlnames.push_back(lname ? lname.Get() : ""); + } + if(newvnames.size() == 0) return "No variables found"; + + { + auto his = newnc.A("history"); + if(his) newhistory = his; + } + + history = std::move(newhistory); + vnames = std::move(newvnames); + lnames = std::move(newlnames); + nc = std::move(newnc); + return ""; +} + +MString TSCDATAData::Info() const +{ + if(!nc) return ""; + + MString out; + + out += MString("Dimensions: ") + nx + " X " + ny + "\n"; + out += MString("Variables: "); + for(size_t i = 0; i < vnames.size(); i++) out += ((i == 0) ? "" : ", ") + vnames[i]; + out += "\n"; + auto his = nc.A("history"); + if(his) out += "Creator: " + his.Get() + "\n"; + return out; +} + +std::vector TSCDATAData::ReadLons() const +{ + std::vector out; + auto lons = nc.V("longitude"); + if(lons) + { + out.resize(lons.DimLen(0)); + for(size_t i = 0; i < out.size(); i++) out[i] = lons(i); + } + return out; +} + +std::vector TSCDATAData::ReadLats() const +{ + std::vector out; + auto lats = nc.V("latitude"); + if(lats) + { + out.resize(lats.DimLen(0)); + for(size_t i = 0; i < out.size(); i++) out[i] = lats(i); + } + return out; +} + +std::vector TSCDATAData::ReadVar(const MString& name) const +{ + std::vector out; + bool havevar = false; + for(size_t i = 0; i < vnames.size(); i++) havevar = havevar || (vnames[i] == name); + if(!havevar) return out; + + auto var = nc.V(name, "latitude", "longitude"); + if(var) + { + out.resize(nx * ny); + for(size_t iy = 0; iy < ny; iy++) + for(size_t ix = 0; ix < nx; ix++) out[iy * nx + ix] = var(iy, ix); + } + return out; +} + +TSCDATAData::DataType TSCDATAData::FillVal(const MString& name) const +{ + auto fill = nc.A(name, "_FillValue"); + return fill ? fill.Get() : 0.0; +} diff --git a/sources/TSCDATA.h b/sources/TSCDATA.h new file mode 100644 index 0000000..080e1b3 --- /dev/null +++ b/sources/TSCDATA.h @@ -0,0 +1,35 @@ +#pragma once +#include "DataAdapters/ncfilealt.h" +#include "ParseArgs.h" +#include "simple2ddata.h" + +class TSCDATAData +{ + michlib::NCFileA nc; + std::vector vnames, lnames; + size_t nx, ny; + MString history; + + using DataType = float; + + public: + static constexpr const char* name = "TSCDATA"; + + MString Info() const; + MString Open(const CLArgs& args); + + std::vector ReadLons() const; + std::vector ReadLats() const; + std::vector ReadVar(const MString& name) const; + + const auto& VarNames() const { return vnames; } + const auto& LongNames() const { return lnames; } + + size_t Nx() const { return nx; } + size_t Ny() const { return ny; } + size_t NVar() const { return vnames.size(); } + + const auto& History() const { return history; } + + DataType FillVal(const MString& name) const; +};