Browse Source

Added support for netcdf output for the tsc action

interpolate
Michael Uleysky 10 months ago
parent
commit
17b04d63de
  1. 158
      actions/actiontsc.h
  2. 2
      include/ncfuncs.h
  3. 14
      include/simple2ddata.h
  4. 10
      include/traits.h
  5. 42
      src/ncfuncs.cpp

158
actions/actiontsc.h

@ -1,6 +1,7 @@
#pragma once
#include "BFileW.h"
#include "actiondep.h"
#include "ncfuncs.h"
#include <memory>
using michlib::BFileW;
@ -49,23 +50,150 @@ template<class D> 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<float>::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<D>)
{
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<D>)
{
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;
};

2
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<MDateTime, time_t, bool> Refdate(const MString& refdate);
static void GetVars(const NCFileA& nc, std::set<MString>& vars);
static CoordNames GetCNames(const NCFileA& nc);

14
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; }
};

10
include/traits.h

@ -130,3 +130,13 @@ concept ReadIsGrid = requires {
ReadType<T>().YStep()
} -> std::convertible_to<real>;
};
template<class T>
concept ReadIs2DGeoRectArray = requires {
{
ReadType<T>().Ix2Lon(0)
} -> std::convertible_to<real>;
{
ReadType<T>().Iy2Lat(0)
} -> std::convertible_to<real>;
};

42
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<class HV> 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;
}
*/

Loading…
Cancel
Save