odm/sources/AVISOLOCAL.cpp

#define MICHLIB_NOSOURCE
#include "AVISOLOCAL.h"

MString AVISOLOCALData::Info() const
{
 if(!isOk()) return "";

 NCFileA           nc;
 struct CoordNames cn, dn;
 std::set<MString> vars;

 nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc");
 if(!nc) return "Can't open file " + datapath + "/uv-" + times[0].ToString() + ".nc";

 GetVars(nc, vars);
 dn = GetDNames(nc);
 cn = GetCNames(nc);

 MString svars;
 {
  bool first = true;
  for(const auto& v: vars)
  {
   svars += (first ? "" : ", ") + v;
   first = false;
  }
 }

 auto lons = nc.VR(cn.lonname);
 auto lats = nc.VR(cn.latname);
 if(!(lons && lats)) return "Can't get longitudes/latitudes";

 real lonb    = lons(0);
 real latb    = lats(0);
 real lone    = lons(dn.nx - 1);
 real late    = lats(dn.ny - 1);
 real lonstep = (lone - lonb) / (dn.nx - 1);
 real latstep = (late - latb) / (dn.ny - 1);

 // clang-format off
 return
  "Dataset: " + DataTitle() + "\n" +
  " Begin date: " + Time(0).ToString() + "\n" +
  " End date: " + Time(NTimes()-1).ToString() + "\n" +
  " Time step: " + Timestep() + " seconds\n" +
  " Time moments: " + NTimes() + "\n" +
  " Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" +
  " Grid: " + dn.nx + "x" + dn.ny + " (" + lonstep + " x " + latstep + ")\n" +
  " Supported variables: " + svars;
 // clang-format on
}

MString AVISOLOCALData::Open(const CLArgs& args)
{
 GPL.UsePrefix("AVISOLOCAL");
 datapath = GPL.ParameterSValue("Datapath", "");
 RegExp regex("uv-([0-9]{4}-[0-9]{2}-[0-9]{2}).nc");
 regex.Compile();

 DIR*           dir = opendir(datapath.Buf());
 struct dirent* de;
 if(nullptr == dir) return "Can't open directory " + datapath;
 while((de = readdir(dir)))
 {
  if(!regex.Match(de->d_name)) continue;
  times.emplace_back(MString(de->d_name + regex.Off(1), regex.Len(1)));
 }
 closedir(dir);
 std::sort(times.begin(), times.end());
 return "";
}

std::pair<const BaseParameters*, MString> AVISOLOCALData::Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const
{
 std::unique_ptr<struct Parameters> ppar{new struct Parameters};
 ppar->lonb = reg.lonb;
 ppar->lone = reg.lone;
 ppar->latb = reg.latb;
 ppar->late = reg.late;

 pars.SetParameter("lonb", ppar->lonb);
 pars.SetParameter("latb", ppar->latb);
 pars.SetParameter("lone", ppar->lone);
 pars.SetParameter("late", ppar->late);

 return {ppar.release(), ""};
}

AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParameters* ip, size_t i) const
{
 if(!isOk()) return Data();

 auto    p = dynamic_cast<const struct Parameters*>(ip);
 NCFileA nc;
 MString name    = "";
 bool    isfloat = false, isint2 = false, isint4 = false;

 nc.Reset(datapath + "/uv-" + times[i].ToString() + ".nc");
 if(!nc) return Data();
 {
  auto head = nc.Header();
  for(const auto& v: head.Variables())
  {
   auto stname = nc.A<MString>(v.Name(), "standard_name");
   if(!stname) continue;
   if(StName2Name(stname) == vname)
   {
    name    = v.Name();
    isint2  = v.Type().Id() == NC_SHORT;
    isint4  = v.Type().Id() == NC_INT;
    isfloat = v.Type().Id() == NC_FLOAT;
   }
  }
 }

 if(!name.Exist()) // Conversion read
 {
  // U and U2 from u and v
  if(vname == "U" || vname == "U2")
  {
   bool square = vname == "U2";
   auto u      = Read("u", ip, i);
   auto v      = Read("v", ip, i);
   if(!(u && v)) return Data();
   auto out = u;
   out.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
   for(size_t ind = 0; ind < out.N(); ind++)
   {
    if(u.IsFill(ind) || v.IsFill(ind))
     out.V(ind) = out.Fillval();
    else
     out.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
   }
   return out;
  }
  return Data();
 }

 // Direct read
 if(isint2) return ReadVarRaw<int2>(nc, name, p);
 if(isint4) return ReadVarRaw<int4>(nc, name, p);
 if(isfloat) return ReadVarRaw<float>(nc, name, p);
 return Data();
}

template<class DataType> AVISOLOCALData::Data AVISOLOCALData::ReadVarRaw(const NCFileA& nc, const MString& name, const struct AVISOLOCALData::Parameters* p) const
{
 real     unitmul = 1.0;
 DataType fill;
 real     offset = 0.0, scale = 1.0;

 {
  auto a_fill     = nc.A<DataType>(name, "_FillValue");
  auto a_offset_d = nc.A<double>(name, "add_offset");
  auto a_scale_d  = nc.A<double>(name, "scale_factor");
  auto a_offset_f = nc.A<float>(name, "add_offset");
  auto a_scale_f  = nc.A<float>(name, "scale_factor");
  if(!a_fill) return Data();
  fill = a_fill;
  if(a_offset_d) offset = a_offset_d;
  if(a_scale_d) scale = a_scale_d;
  if(a_offset_f) offset = a_offset_f;
  if(a_scale_f) scale = a_scale_f;
 }

 auto cn = GetCNames(nc);
 auto dn = GetDNames(nc);

 auto lons = nc.VR(cn.lonname);
 auto lats = nc.VR(cn.latname);
 if(!(lons && lats)) return Data();

 auto lonb    = lons(0);
 auto latb    = lats(0);
 auto lone    = lons(dn.nx - 1);
 auto late    = lats(dn.ny - 1);
 auto lonstep = (lone - lonb) / (dn.nx - 1);
 auto latstep = (late - latb) / (dn.ny - 1);

 auto dom  = DetGeoDomain(lonb, lone);
 real lon1 = ToGeoDomain(p->lonb, dom);
 real lon2 = ToGeoDomain(p->lone, dom);
 real lat1 = p->latb;
 real lat2 = p->late;

 // Special case when the longitude lies in a small sector between the end and the start
 if(lon1 < lonb) lon1 = lone;
 if(lon2 > lone) lon2 = lonb;

 auto yb = static_cast<size_t>(Floor((lat1 - latb) / latstep));
 auto ye = static_cast<size_t>(Ceil((lat2 - latb) / latstep));
 if(ye > dn.ny - 1) ye = dn.ny - 1;
 if(yb >= ye) return Data();

 auto xb = static_cast<size_t>(Floor((lon1 - lonb) / lonstep));
 auto xe = static_cast<size_t>(Ceil((lon2 - lonb) / lonstep));

 if(xb == xe) return Data();

 MString unit;
 {
  auto unitatt = nc.A<MString>(name, "units");
  if(unitatt) unit = unitatt;
 }
 if(unit == "m s-1" || unit == "m/s")
 {
  unitmul = 100.0;
  unit    = "cm/s";
 }

 Data data((xb < xe) ? (xe - xb + 1) : (dn.nx + xe - xb + 1), ye - yb + 1, lons(xb), lats(yb), lonstep, latstep, std::move(unit));

 if(xb < xe)
 {
  auto var = nc.V<DataType>(name, {dn.lonname, xb, xe - xb + 1}, {dn.latname, yb, ye - yb + 1});

  if(!var) return Data();
  if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data();

  for(size_t ix = 0; ix < var.DimLen(0); ix++)
   for(size_t iy = 0; iy < var.DimLen(1); iy++)
   {
    DataType v   = var(ix, iy);
    data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
   }
 }
 else
 {
  auto var1 = nc.V<DataType>(name, {dn.lonname, xb}, {dn.latname, yb, ye - yb + 1});
  auto var2 = nc.V<DataType>(name, {dn.lonname, 0, xe + 1}, {dn.latname, yb, ye - yb + 1});

  if(!(var1 && var2)) return Data();
  if((var1.DimLen(0) + var2.DimLen(0)) != data.Nx() || var1.DimLen(1) != data.Ny() || var2.DimLen(1) != data.Ny()) return Data();
  for(size_t ix = 0; ix < var1.DimLen(0); ix++)
   for(size_t iy = 0; iy < var1.DimLen(1); iy++)
   {
    DataType v   = var1(ix, iy);
    data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
   }
  for(size_t ix = 0; ix < var2.DimLen(0); ix++)
   for(size_t iy = 0; iy < var2.DimLen(1); iy++)
   {
    DataType v                    = var2(ix, iy);
    data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
   }
 }
 return data;
}