odm/actions/actiongrad.cpp

#define MICHLIB_NOSOURCE
#include "actiongrad.h"

MString GradMethods::NCFileW::Create(const MString& name, const MString& history, const std::vector<MString>& vnames, const std::vector<MString>& lnames,
                                     const std::vector<GradMethods::DataType>& lons, const std::vector<GradMethods::DataType>& 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<GradMethods::DataType>(std::numeric_limits<GradMethods::Matrix::MDataType>::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<GradMethods::Matrix::MDataType>& 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<GradMethods::DataType>(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<GradMethods::DataType>& 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<sizeof(DataType)>::Log(minmax.min);
  minmax.max = michlib_internal::RealType<sizeof(DataType)>::Log(minmax.max);
 }

 DataType a = (std::numeric_limits<MDataType>::max() - 1) / (minmax.max - minmax.min);
 for(size_t i = 1; i < in.size(); i++)
 {
  DataType v = minmax.log ? michlib_internal::RealType<sizeof(DataType)>::Log(in[i]) : in[i];
  if(in[i] == minmax.fill)
   data[i] = std::numeric_limits<MDataType>::max();
  else if(v <= minmax.min)
   data[i] = 0;
  else if(v >= minmax.max)
   data[i] = std::numeric_limits<MDataType>::max() - 1;
  else
   data[i] = static_cast<MDataType>(michlib_internal::RealType<sizeof(DataType)>::Round(a * (v - minmax.min)));
 }
}

void GradMethods::Matrix::Grad()
{
 std::vector<MDataType> out(data.size());
 const auto             bad = std::numeric_limits<MDataType>::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<IT>(michlib::Round(michlib::Hypot(gx, gy)));
    if(sq >= bad) sq = bad - 1;
    out[iy * nx + ix] = static_cast<MDataType>(sq);
   }
  }

 data = std::move(out);
}