diff --git a/actions/actiontsc.h b/actions/actiontsc.h index 6726a97..37d1ffa 100644 --- a/actions/actiontsc.h +++ b/actions/actiontsc.h @@ -1,6 +1,7 @@ #pragma once #include "BFileW.h" #include "actiondep.h" +#include "ncfuncs.h" #include using michlib::BFileW; @@ -49,23 +50,150 @@ template MString ActionTSC::DoAction(const CLArgs& args, D& ds) 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") : ""; - if(!name.Exist()) name = "out.bin"; + MString outfmt = args.contains("format") ? args.at("format") : "bin"; - fw.Create(name, 3); - fw.SetColumnName(1, "Longitude"); - fw.SetColumnName(2, "Latitude"); - fw.SetColumnName(3, vname + ", " + (data.Unit().Exist() ? data.Unit() : "unknown")); - 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(i) == data.Fillval() ? NAN : data(i)); + BFileW fw; + MString name = args.contains("out") ? args.at("out") : "out.bin"; + + fw.Create(name, 3); + fw.SetColumnName(1, "Longitude"); + fw.SetColumnName(2, "Latitude"); + fw.SetColumnName(3, vname + ", " + (data.Unit().Exist() ? data.Unit() : "unknown")); + 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.IsFill(i) ? NAN : data(i)); + } + fw.Finalize(); + fw.Close(); + return ""; + } + + if(outfmt == "nc" || outfmt == "netcdf") + { + MString name = args.contains("out") ? args.at("out") : "out.nc"; + int ret; + int ncid, dimid, dimidxy[2]; + int lonid, latid, vid; + MString text; + auto fill = std::numeric_limits::max(); + float fval; + int compress = 3; + + if(args.contains("compress")) compress = args.at("compress").ToInt(); + + 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", args.at("_cmdline").Len() + 1, args.at("_cmdline").Buf()); + if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file"; + + if constexpr(ReadIs2DGeoRectArray) + { + ret = nc_def_dim(ncid, "longitude", data.Nx(), &dimidxy[0]); + if(ret != NC_NOERR) return "Can't create x-dimension in the netcdf file"; + ret = nc_def_dim(ncid, "latitude", data.Ny(), &dimidxy[1]); + if(ret != NC_NOERR) return "Can't create y-dimension in the netcdf file"; + ret = nc_def_var(ncid, "longitude", NC_FLOAT, 1, &dimidxy[0], &lonid); + if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file"; + ret = nc_def_var(ncid, "latitude", NC_FLOAT, 1, &dimidxy[1], &latid); + if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file"; + ret = nc_def_var(ncid, vname.Buf(), NC_FLOAT, 2, dimidxy, &vid); + if(ret != NC_NOERR) return "Can't create " + vname + " variable in the netcdf file"; + } + else + { + ret = nc_def_dim(ncid, "i", data.N(), &dimid); + if(ret != NC_NOERR) return "Can't create dimension in the netcdf file"; + ret = nc_def_var(ncid, "longitude", NC_FLOAT, 1, &dimid, &lonid); + if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file"; + ret = nc_def_var(ncid, "latitude", NC_FLOAT, 1, &dimid, &latid); + if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file"; + ret = nc_def_var(ncid, vname.Buf(), NC_FLOAT, 1, &dimid, &vid); + if(ret != NC_NOERR) return "Can't create " + vname + " variable in the netcdf file"; + } + + ret = nc_def_var_deflate(ncid, lonid, 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(ncid, latid, 1, 1, compress); + if(ret != NC_NOERR) return "Can't set deflate parameters for latitude variable in the netcdf file"; + ret = nc_def_var_deflate(ncid, vid, 1, 1, compress); + if(ret != NC_NOERR) return "Can't set deflate parameters for " + vname + " variable in the netcdf file"; + + text = "longitude"; + ret = nc_put_att_text(ncid, lonid, "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, latid, "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 = NCFuncs::Name2StName(vname); + if(text.Exist()) + { + ret = nc_put_att_text(ncid, vid, "standard_name", text.Len() + 1, text.Buf()); + if(ret != NC_NOERR) return "Can't write standard_name attribute of " + vname + " variable in the netcdf file"; + } + else + { + text = NCFuncs::Name2LongName(vname); + ret = nc_put_att_text(ncid, vid, "long_name", text.Len() + 1, text.Buf()); + if(ret != NC_NOERR) return "Can't write long_name attribute of " + vname + " variable in the netcdf file"; + } + + ret = nc_put_att_float(ncid, vid, "_FillValue", NC_FLOAT, 1, &fill); + if(ret != NC_NOERR) return "Can't write _FillValue attribute of " + vname + " variable in the netcdf file"; + + ret = nc_enddef(ncid); + if(ret != NC_NOERR) return "Can't finish definition of the netcdf file"; + + if constexpr(ReadIs2DGeoRectArray) + { + size_t i[2]; + for(i[0] = 0; i[0] < data.Nx(); i[0]++) + { + fval = data.Ix2Lon(i[0]); + ret = nc_put_var1_float(ncid, lonid, &i[0], &fval); + if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file: " + MString(i[0]); + } + for(i[1] = 0; i[1] < data.Ny(); i[1]++) + { + fval = data.Iy2Lat(i[1]); + ret = nc_put_var1_float(ncid, latid, &i[1], &fval); + if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file: " + MString(i[1]); + } + for(i[0] = 0; i[0] < data.Nx(); i[0]++) + for(i[1] = 0; i[1] < data.Ny(); i[1]++) + { + fval = data.IsFill(i[0], i[1]) ? fill : data(i[0], i[1]); + ret = nc_put_var1_float(ncid, vid, i, &fval); + if(ret != NC_NOERR) return "Can't write " + vname + " variable in the netcdf file: " + MString(i[0]) + ", " + i[1]; + } + } + else + { + for(size_t i = 0; i < data.N(); i++) + { + fval = data.Lon(i); + ret = nc_put_var1_float(ncid, lonid, &i, &fval); + if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file: " + MString(i); + fval = data.Lat(i); + ret = nc_put_var1_float(ncid, latid, &i, &fval); + if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file: " + MString(i); + fval = data.IsFill(i) ? fill : data(i); + ret = nc_put_var1_float(ncid, vid, &i, &fval); + if(ret != NC_NOERR) return "Can't write " + vname + " variable in the netcdf file: " + MString(i); + } + } + + ret = nc_close(ncid); + if(ret != NC_NOERR) return "Can't close netcdf file"; + + return ""; } - fw.Finalize(); - fw.Close(); - return ""; + return "Unknown format: " + outfmt; }; diff --git a/include/ncfuncs.h b/include/ncfuncs.h index f77be7b..2dcf444 100644 --- a/include/ncfuncs.h +++ b/include/ncfuncs.h @@ -18,6 +18,8 @@ class NCFuncs }; static MString StName2Name(const MString& stname); + static MString Name2StName(const MString& name); + static MString Name2LongName(const MString& name); static std::tuple Refdate(const MString& refdate); static void GetVars(const NCFileA& nc, std::set& vars); static CoordNames GetCNames(const NCFileA& nc); diff --git a/include/simple2ddata.h b/include/simple2ddata.h index 80bc958..66b9924 100644 --- a/include/simple2ddata.h +++ b/include/simple2ddata.h @@ -8,6 +8,9 @@ class Simple2DData: public BaseData real xstep = 0.0, ystep = 0.0; public: + using BaseData::IsFill; + using BaseData::V; + using BaseData::operator(); Simple2DData() = default; Simple2DData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_, MString&& unit = ""): @@ -15,12 +18,6 @@ class Simple2DData: public BaseData { } - 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); } @@ -30,12 +27,17 @@ class Simple2DData: public BaseData size_t Nx() const { return nx; } size_t Ny() const { return ny; } + bool IsFill(size_t ix, size_t iy) const { return IsFill(iy * nx + ix); } + real Lon(size_t ix, [[maybe_unused]] size_t iy) const { return x0 + ix * xstep; } real Lat([[maybe_unused]] size_t ix, size_t iy) const { return y0 + iy * ystep; } real Lon(size_t i) const { return Lon(i % nx, i / nx); } real Lat(size_t i) const { return Lat(i % nx, i / nx); } + real Ix2Lon(size_t ix) const { return x0 + ix * xstep; } + real Iy2Lat(size_t iy) const { return y0 + iy * ystep; } + real XStep() const { return xstep; } real YStep() const { return ystep; } }; diff --git a/include/traits.h b/include/traits.h index 29271fa..b0d85d7 100644 --- a/include/traits.h +++ b/include/traits.h @@ -130,3 +130,13 @@ concept ReadIsGrid = requires { ReadType().YStep() } -> std::convertible_to; }; + +template +concept ReadIs2DGeoRectArray = requires { + { + ReadType().Ix2Lon(0) + } -> std::convertible_to; + { + ReadType().Iy2Lat(0) + } -> std::convertible_to; +}; diff --git a/src/ncfuncs.cpp b/src/ncfuncs.cpp index 054ce3f..7d514c9 100644 --- a/src/ncfuncs.cpp +++ b/src/ncfuncs.cpp @@ -126,6 +126,32 @@ MString NCFuncs::StName2Name(const MString& stname) return ""; } +MString NCFuncs::Name2StName(const MString& name) +{ + if(name == "ptemp") return "sea_water_potential_temperature"; + if(name == "temp") return "sea_water_temperature"; + if(name == "sal") return "sea_water_salinity"; + if(name == "mld") return "ocean_mixed_layer_thickness_defined_by_sigma_theta"; + if(name == "ssh") return "sea_surface_height_above_geoid"; + if(name == "u") return "eastward_sea_water_velocity"; + if(name == "v") return "northward_sea_water_velocity"; + if(name == "w") return "upward_sea_water_velocity"; + if(name == "chl") return "mass_concentration_of_chlorophyll_a_in_sea_water"; + if(name == "NO3") return "mole_concentration_of_nitrate_in_sea_water"; + if(name == "prprod") return "net_primary_production_of_biomass_expressed_as_carbon_per_unit_volume_in_sea_water"; + if(name == "Cchl") return "mole_concentration_of_phytoplankton_expressed_as_carbon_in_sea_water"; + if(name == "PO4") return "mole_concentration_of_phosphate_in_sea_water"; + if(name == "Si") return "mole_concentration_of_silicate_in_sea_water"; + return ""; +} + +MString NCFuncs::Name2LongName(const MString& name) +{ + if(name == "U") return "Module of horizontal velocity"; + if(name == "U2") return "Squared horizontal velocity"; + return ""; +} + bool NCFuncs::HaveVar(const NCFileA& nc, const MString& vname) { auto head = nc.Header(); @@ -137,19 +163,3 @@ bool NCFuncs::HaveVar(const NCFileA& nc, const MString& vname) } return false; } - -/* -template bool NCFuncs::CheckVar(const MString& vname, HV hv) -{ - if(!hv(vname)) - { - bool varexist = false; - if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = true; - if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = true; - if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = true; - if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = true; - if(!varexist) return false; - } - return true; -} -*/