From ed1282413f52afad7d214fe342b2adefad7af904 Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Sun, 30 Jul 2023 17:29:50 +1000 Subject: [PATCH] Added support for netcdf output for the uv action --- actions/actionuv.cpp | 88 ++++++++++++++ actions/actionuv.h | 277 +++++++++++++++++++++++++++++++++---------- 2 files changed, 302 insertions(+), 63 deletions(-) create mode 100644 actions/actionuv.cpp diff --git a/actions/actionuv.cpp b/actions/actionuv.cpp new file mode 100644 index 0000000..434d72c --- /dev/null +++ b/actions/actionuv.cpp @@ -0,0 +1,88 @@ +#define MICHLIB_NOSOURCE +#include "actionuv.h" + +void UVMethods::StPoints::WriteBinBile(const MString& name) const +{ + BFileW stp; + stp.Create(name, 3); + stp.SetColumnName(1, lonlat ? "Longitude" : "x"); + stp.SetColumnName(2, lonlat ? "Latitude" : "x"); + stp.SetColumnName(3, "Stability (0 - saddle, 1 - st. anticicl. focus, 2 - st. knot, 3 - unst. anticicl. focus, 4 - unst. knot, 5 - st. cicl. focus, 6 - unst. cicl. focus)"); + + for(size_t i = 0; i < N(); i++) + { + stp.Write(x[i]); + stp.Write(y[i]); + stp.Write(michlib::int_cast(t[i])); + } + stp.Finalize(); + stp.Close(); +} + +MString UVMethods::StPoints::WriteNcFile(const MString& name, const MString& history, int comp) const +{ + int ret; + int ncid; + int dimid; + int xid, yid, tid; + MString text; + + ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &ncid); + if(ret != NC_NOERR) return "Can't create netcdf file: " + name; + + ret = nc_put_att_text(ncid, NC_GLOBAL, "history", history.Len() + 1, history.Buf()); + if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file"; + + ret = nc_def_dim(ncid, "i", N(), &dimid); + if(ret != NC_NOERR) return "Can't create dimension in the netcdf file"; + + ret = nc_def_var(ncid, lonlat ? "longitude" : "x", NC_FLOAT, 1, &dimid, &xid); + if(ret != NC_NOERR) return "Can't create x variable in the netcdf file"; + ret = nc_def_var(ncid, lonlat ? "latitude" : "y", NC_FLOAT, 1, &dimid, &yid); + if(ret != NC_NOERR) return "Can't create y variable in the netcdf file"; + + if(lonlat) + { + text = "longitude"; + ret = nc_put_att_text(ncid, 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(ncid, 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"; + } + + ret = nc_def_var_deflate(ncid, xid, 1, 1, comp); + if(ret != NC_NOERR) return "Can't set deflate parameters for x variable in the netcdf file"; + ret = nc_def_var_deflate(ncid, yid, 1, 1, comp); + if(ret != NC_NOERR) return "Can't set deflate parameters for y variable in the netcdf file"; + + ret = nc_def_var(ncid, "type", NC_UBYTE, 1, &dimid, &tid); + if(ret != NC_NOERR) return "Can't create type variable in the netcdf file"; + ret = nc_def_var_deflate(ncid, tid, 1, 1, comp); + if(ret != NC_NOERR) return "Can't set deflate parameters for type variable in the netcdf file"; + text = "Stationary point type, 0 - saddle, 1 - st. anticicl. focus, 2 - st. knot, 3 - unst. anticicl. focus, 4 - unst. knot, 5 - st. cicl. focus, 6 - unst. cicl. focus"; + ret = nc_put_att_text(ncid, tid, "long_name", text.Len() + 1, text.Buf()); + if(ret != NC_NOERR) return "Can't write long_name attribute of type variable in the netcdf file"; + + ret = nc_enddef(ncid); + if(ret != NC_NOERR) return "Can't finish definition of the netcdf file"; + + const size_t i = 0, c = N(); + float buf[c]; + + for(size_t ix = 0; ix < c; ix++) buf[ix] = x[ix]; + ret = nc_put_vara_float(ncid, xid, &i, &c, buf); + if(ret != NC_NOERR) return "Can't write x variable in the netcdf file"; + + for(size_t ix = 0; ix < c; ix++) buf[ix] = y[ix]; + ret = nc_put_vara_float(ncid, yid, &i, &c, buf); + if(ret != NC_NOERR) return "Can't write y variable in the netcdf file"; + + ret = nc_put_vara_ubyte(ncid, tid, &i, &c, t.data()); + if(ret != NC_NOERR) return "Can't write type variable in the netcdf file"; + + ret = nc_close(ncid); + if(ret != NC_NOERR) return "Can't close netcdf file"; + + return ""; +} diff --git a/actions/actionuv.h b/actions/actionuv.h index da64ebf..810b267 100644 --- a/actions/actionuv.h +++ b/actions/actionuv.h @@ -1,11 +1,98 @@ #pragma once #include "BFileW.h" #include "actiondep.h" -#include +#include "ncfilew.h" +#include using michlib::BFileW; -ADD_ACTION(UV, uv, ReadPSupported || ReadSupported); +class UVMethods +{ + protected: + template class SparserHolder + { + protected: + const D& data; + size_t shiftx, shifty, skipx, skipy; + SparserHolder(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): data(d), shiftx(shx), shifty(shy), skipx(sx), skipy(sy) {} + + size_t I2Ix(size_t i) const { return i * skipx + shiftx; } + size_t I2Iy(size_t i) const { return i * skipy + shifty; } + + public: + size_t Nx() const + { + size_t n = data.Nx() - shiftx; + return n / skipx + ((n % skipx > 0) ? 1 : 0); + } + size_t Ny() const + { + size_t n = data.Ny() - shifty; + return n / skipy + ((n % skipy > 0) ? 1 : 0); + } + + static constexpr auto Fillval() { return D::Fillval(); } + + real Lon(size_t ix, size_t iy) const { return data.Lon(I2Ix(ix), I2Iy(iy)); } + real Lat(size_t ix, size_t iy) const { return data.Lat(I2Ix(ix), I2Iy(iy)); } + real U(size_t ix, size_t iy) const { return data.U(I2Ix(ix), I2Iy(iy)); } + real V(size_t ix, size_t iy) const { return data.V(I2Ix(ix), I2Iy(iy)); } + bool IsFill(size_t ix, size_t iy) const { return data.IsFill(I2Ix(ix), I2Iy(iy)); } + }; + + template class Sparser2DGeoRectArray: public SparserHolder + { + protected: + Sparser2DGeoRectArray(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): SparserHolder(d, shx, shy, sx, sy) {} + }; + template class Sparser2DGeoRectArray: public SparserHolder + { + protected: + using SparserHolder::data; + using SparserHolder::I2Ix; + using SparserHolder::I2Iy; + + Sparser2DGeoRectArray(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): SparserHolder(d, shx, shy, sx, sy) {} + + public: + real Ix2Lon(size_t ix) const { return data.Ix2Lon(I2Ix(ix)); } + real Iy2Lat(size_t iy) const { return data.Iy2Lat(I2Iy(iy)); } + }; + + template class Sparser: public Sparser2DGeoRectArray> + { + public: + Sparser(const D& d, size_t shx, size_t shy, size_t sx, size_t sy): Sparser2DGeoRectArray>(d, shx, shy, sx, sy) {} + }; + + class StPoints + { + std::vector x, y; + std::vector t; + bool lonlat; + + size_t N() const { return t.size(); } + + public: + StPoints(bool ll): lonlat(ll) {} + + void Add(const std::vector& points) + { + for(const auto& p: points) + { + x.push_back(p.x); + y.push_back(p.y); + t.push_back(michlib::int_cast(p.type)); + } + } + + void WriteBinBile(const MString& name) const; + + MString WriteNcFile(const MString& name, const MString& history, int comp) const; + }; +}; + +ADD_ACTION(UV, uv, ReadPSupported || ReadSupported, UVMethods); template MString ActionUV::DoAction(const CLArgs& args, D& ds) { @@ -44,39 +131,81 @@ template MString ActionUV::DoAction(const CLArgs& args, D& ds) if(!data) return "Can't read data"; // Main file - MString name = args.contains("out") ? args.at("out") : ""; + MString name = args.contains("out") ? args.at("out") : ""; + MString outfmt = args.contains("outformat") ? args.at("outformat") : (GetExt(name) == "nc" ? "nc" : "bin"); + + MString u = data.Unit().Exist() ? data.Unit() : "unknown", d = data.DUnit().Exist() ? data.DUnit() : "unknown"; + if(name.Exist()) { - BFileW fw; - fw.Create(name, 9); - fw.SetColumnName(1, "Longitude"); - fw.SetColumnName(2, "Latitude"); - fw.SetColumnName(3, "u, cm/s"); - fw.SetColumnName(4, "v, cm/s"); - fw.SetColumnName(5, "Div, (cm/s)/km"); - fw.SetColumnName(6, "Rot, (cm/s)/km"); - fw.SetColumnName(7, "Okubo-Weiss parameter, (cm^2/s^2)/km^2"); - fw.SetColumnName(8, "Kinetic energy, cm^2/s^2"); - fw.SetColumnName(9, "Eddy kinetic energy, cm^2/s^2"); - fw.SetParameters(pars); - for(size_t i = 0; i < data.N(); i++) + if(outfmt == "bin") { - fw.Write(data.Lon(i)); - fw.Write(data.Lat(i)); - fw.Write(data.U(i) == data.Fillval() ? NAN : data.U(i)); - fw.Write(data.V(i) == data.Fillval() ? NAN : data.V(i)); - fw.Write(data.Div(i) == data.Fillval() ? NAN : data.Div(i)); - fw.Write(data.Rot(i) == data.Fillval() ? NAN : data.Rot(i)); - fw.Write(data.OW(i) == data.Fillval() ? NAN : data.OW(i)); - fw.Write(data.U2(i) == data.Fillval() ? NAN : data.U2(i)); - fw.Write((data.U(i) == data.Fillval() || data.V(i) == data.Fillval()) ? NAN : (data.U2(i) - (data.U(i) * data.U(i) + data.V(i) * data.V(i)))); + BFileW fw; + fw.Create(name, 9); + fw.SetColumnName(1, "Longitude"); + fw.SetColumnName(2, "Latitude"); + fw.SetColumnName(3, "u, " + u); + fw.SetColumnName(4, "v, " + u); + fw.SetColumnName(5, "Div, (" + u + ")/" + d); + fw.SetColumnName(6, "Rot, (" + u + ")/" + d); + fw.SetColumnName(7, "Okubo-Weiss parameter, (" + u + ")2/" + d + "2"); + fw.SetColumnName(8, "Kinetic energy, (" + u + ")2"); + fw.SetColumnName(9, "Eddy kinetic energy, (" + u + ")2"); + fw.SetParameters(pars); + + for(size_t i = 0; i < data.N(); i++) + { + fw.Write(data.Lon(i)); + fw.Write(data.Lat(i)); + fw.Write(data.U(i) == data.Fillval() ? NAN : data.U(i)); + fw.Write(data.V(i) == data.Fillval() ? NAN : data.V(i)); + fw.Write(data.Div(i) == data.Fillval() ? NAN : data.Div(i)); + fw.Write(data.Rot(i) == data.Fillval() ? NAN : data.Rot(i)); + fw.Write(data.OW(i) == data.Fillval() ? NAN : data.OW(i)); + fw.Write(data.U2(i) == data.Fillval() ? NAN : data.U2(i)); + fw.Write((data.U(i) == data.Fillval() || data.V(i) == data.Fillval()) ? NAN : (data.U2(i) - (data.U(i) * data.U(i) + data.V(i) * data.V(i)))); + } + fw.Finalize(); + fw.Close(); } - fw.Finalize(); - fw.Close(); + else if(outfmt == "nc" || outfmt == "netcdf") + { + int compress = 3; + NCFileW fw; + MString err; + + if(args.contains("compress")) compress = args.at("compress").ToInt(); + + if(!err.Exist()) err = fw.Create(data, name, args.at("_cmdline"), compress); + if(!err.Exist()) err = fw.AddVariable("u", "", "Eastward velocity", u); + if(!err.Exist()) err = fw.AddVariable("v", "", "Northward velocity", u); + if(!err.Exist()) err = fw.AddVariable("div", "", "Velocity divergence", "(" + u + ")/" + d); + if(!err.Exist()) err = fw.AddVariable("rot", "", "Velocity rotor", "(" + u + ")/" + d); + if(!err.Exist()) err = fw.AddVariable("ow", "", "Okubo-Weiss parameter", "(" + u + ")2/" + d + "2"); + if(!err.Exist()) err = fw.AddVariable("ke", "", "Squared velocity module, u^2+v^2", "(" + u + ")2"); + if(!err.Exist()) err = fw.AddVariable("eke", "", "Squared velocity dispersion aka eddy kinetic energy, -^2-^2", "(" + u + ")2"); + + if(!err.Exist()) err = fw.WriteGrid(data); + + if(!err.Exist()) err = fw.WriteVariable(data, "u", [&data = std::as_const(data)](size_t i, size_t j) { return data.U(i, j); }); + if(!err.Exist()) err = fw.WriteVariable(data, "v", [&data = std::as_const(data)](size_t i, size_t j) { return data.V(i, j); }); + if(!err.Exist()) err = fw.WriteVariable(data, "div", [&data = std::as_const(data)](size_t i, size_t j) { return data.Div(i, j); }); + if(!err.Exist()) err = fw.WriteVariable(data, "rot", [&data = std::as_const(data)](size_t i, size_t j) { return data.Rot(i, j); }); + if(!err.Exist()) err = fw.WriteVariable(data, "ow", [&data = std::as_const(data)](size_t i, size_t j) { return data.OW(i, j); }); + if(!err.Exist()) err = fw.WriteVariable(data, "ke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j); }); + if(!err.Exist()) + err = fw.WriteVariable(data, "eke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j) - (data.U(i, j) * data.U(i, j) + data.V(i, j) * data.V(i, j)); }); + if(err.Exist()) return err; + + fw.Close(); + } + else + return "Unknown format: " + outfmt; } // Filtered vectors file - name = args.contains("velout") ? args.at("velout") : ""; + name = args.contains("velout") ? args.at("velout") : ""; + outfmt = args.contains("veloutformat") ? args.at("veloutformat") : (GetExt(name) == "nc" ? "nc" : "bin"); if(name.Exist()) { size_t shiftx = args.contains("shiftx") ? args.at("shiftx").ToInteger() : 0; @@ -84,48 +213,70 @@ template MString ActionUV::DoAction(const CLArgs& args, D& ds) size_t skipx = args.contains("skipx") ? args.at("skipx").ToInteger() : 1; size_t skipy = args.contains("skipy") ? args.at("skipy").ToInteger() : 1; - BFileW vel; - vel.Create(name, 4); - vel.SetColumnName(1, "Longitude"); - vel.SetColumnName(2, "Latitude"); - vel.SetColumnName(3, "u, cm/s"); - vel.SetColumnName(4, "v, cm/s"); + Sparser sdata(data, shiftx, shifty, skipx, skipy); + if(outfmt == "bin") + { + BFileW vel; + vel.Create(name, 4); + vel.SetColumnName(1, "Longitude"); + vel.SetColumnName(2, "Latitude"); + vel.SetColumnName(3, "u, " + u); + vel.SetColumnName(4, "v, " + u); - for(size_t ix = shiftx; ix < data.Nx(); ix += skipx) - for(size_t iy = shifty; iy < data.Ny(); iy += skipy) - { - vel.Write(data.Lon(ix, iy)); - vel.Write(data.Lat(ix, iy)); - vel.Write(data.U(ix, iy) == data.Fillval() ? NAN : data.U(ix, iy)); - vel.Write(data.V(ix, iy) == data.Fillval() ? NAN : data.V(ix, iy)); - } - vel.Finalize(); - vel.Close(); + for(size_t ix = 0; ix < sdata.Nx(); ix++) + for(size_t iy = 0; iy < sdata.Ny(); iy++) + { + vel.Write(sdata.Lon(ix, iy)); + vel.Write(sdata.Lat(ix, iy)); + vel.Write(sdata.U(ix, iy) == sdata.Fillval() ? NAN : sdata.U(ix, iy)); + vel.Write(sdata.V(ix, iy) == sdata.Fillval() ? NAN : sdata.V(ix, iy)); + } + vel.Finalize(); + vel.Close(); + } + else if(outfmt == "nc" || outfmt == "netcdf") + { + int compress = 3; + NCFileW fw; + MString err; + + if(args.contains("compress")) compress = args.at("compress").ToInt(); + + if(!err.Exist()) err = fw.Create(sdata, name, args.at("_cmdline"), compress); + if(!err.Exist()) err = fw.AddVariable("u", "", "Eastward velocity", u); + if(!err.Exist()) err = fw.AddVariable("v", "", "Northward velocity", u); + + if(!err.Exist()) err = fw.WriteGrid(sdata); + + if(!err.Exist()) err = fw.WriteVariable(sdata, "u", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.U(i, j); }); + if(!err.Exist()) err = fw.WriteVariable(sdata, "v", [&data = std::as_const(sdata)](size_t i, size_t j) { return data.V(i, j); }); + if(err.Exist()) return err; + + fw.Close(); + } + else + return "Unknown format: " + outfmt; } // Stationary points - name = args.contains("stpout") ? args.at("stpout") : ""; + name = args.contains("stpout") ? args.at("stpout") : ""; + outfmt = args.contains("stpoutformat") ? args.at("stpoutformat") : (GetExt(name) == "nc" ? "nc" : "bin"); if(name.Exist()) { - BFileW stp; - stp.Create(name, 3); - stp.SetColumnName(1, "Longitude"); - stp.SetColumnName(2, "Latitude"); - stp.SetColumnName(3, "Stability (0 - saddle, 1 - st. anticicl. focus, 2 - st. knot, 3 - unst. anticicl. focus, 4 - unst. knot, 5 - st. cicl. focus, 6 - unst. cicl. focus)"); - + StPoints p(true); for(size_t ix = 0; ix < data.Nx() - 1; ix++) - for(size_t iy = 0; iy < data.Ny() - 1; iy++) - { - auto points = data.StablePoints(ix, iy); - for(size_t i = 0; i < points.size(); i++) - { - stp.Write(points[i].x); - stp.Write(points[i].y); - stp.Write(michlib::int_cast(points[i].type)); - } - } - stp.Finalize(); - stp.Close(); + for(size_t iy = 0; iy < data.Ny() - 1; iy++) p.Add(data.StablePoints(ix, iy)); + + if(outfmt == "bin") + p.WriteBinBile(name); + else if(outfmt == "nc" || outfmt == "netcdf") + { + int compress = args.contains("compress") ? args.at("compress").ToInt() : 3; + MString err = p.WriteNcFile(name, args.at("_cmdline"), compress); + if(err.Exist()) return err; + } + else + return "Unknown format: " + outfmt; } return "";