odm/src/layereddata.cpp

#define MICHLIB_NOSOURCE
#include "layereddata.h"

MString LayeredData::Info() const
{
 if(!isOk()) return "";
 MString d;
 for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")";

 std::set<MString> vars;
 for(const auto& f: nc)
 {
  auto head = f->Header();
  for(const auto& v: head.Variables())
  {
   auto ret = f->A<MString>(v.Name(), "standard_name");
   if(!ret) continue;
   if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret));
  }
 }
 if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens");
 if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp");
 if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp");

 if(vars.contains("u") && vars.contains("v")) vars.emplace("U");
 if(vars.contains("u") && vars.contains("v")) vars.emplace("U2");

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

 // clang-format off
 return
  "Dataset: " + Title() + "\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: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" +
  " Depths:" + d + "\n" +
  " Supported variables: " + svars;
 // clang-format on
}

MString LayeredData::Open(const MString& dataset)
{
 nc.clear();
 MString proxyurl = GPL.ParameterSValue("USEPROXY", "");
 if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl);

 nc.clear();
 size_t i = 1;
 while(true)
 {
  MString url = GPL.ParameterSValue(dataset + "_URL" + i, "");
  if(url.Exist())
  {
   //michlib::message("Open "+url);
   nc.emplace_back(std::move(url));
   if(!nc.back())
   {
    nc.clear();
    return "Can't connect to url " + url;
   }
   MString ret = nc.back().ReadTimes();
   if(ret.Exist())
   {
    nc.clear();
    return "Can't connect to url " + url + ": " + ret;
   }
  }
  else
   break;
  i++;
 }
 if(nc.size() == 0) return "No urls for dataset " + dataset + " specified in config";

 auto head = nc[0]->Header();
 lonname = latname = "";
 for(const auto& dim: head.Dimensions())
 {
  if(dim.Name() == "lon" || dim.Name() == "longitude")
  {
   lonname = dim.Name();
   nx      = dim.Len();
  }
  if(dim.Name() == "lat" || dim.Name() == "latitude")
  {
   latname = dim.Name();
   ny      = dim.Len();
  }
 }
 if(!(lonname.Exist() && latname.Exist()))
 {
  nc.clear();
  return "Can't find longitude/latitude";
 }

 // Read times
 {
  for(const auto& f: nc) times.insert(times.end(), f.Times().begin(), f.Times().end());
  std::sort(times.begin(), times.end());
  auto last = std::unique(times.begin(), times.end());
  times.erase(last, times.end());
 }

 auto rdepths = nc[0]->VR("depth");
 if(!rdepths)
 {
  nc.clear();
  return "Can't read depths";
 }
 depths.resize(rdepths.DimLen(0));
 for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i);

 auto lons = nc[0]->VR(lonname);
 auto lats = nc[0]->VR(latname);
 if(!(lons && lats))
 {
  nc.clear();
  return "Can't get longitudes/latitudes";
 }
 lonb    = lons(0);
 latb    = lats(0);
 lone    = lons(nx - 1);
 late    = lats(ny - 1);
 lonstep = (lone - lonb) / (nx - 1);
 latstep = (late - latb) / (ny - 1);

 return "";
}

std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const
{
 std::unique_ptr<struct Parameters> ppar{new struct Parameters};

 if(!args.contains("var")) return {nullptr, "Variable not specified"};
 ppar->varname = args.at("var");
 if(!HaveVar(ppar->varname))
 {
  bool varexist = false;
  if(ppar->varname == "temp" && HaveVar("ptemp") && HaveVar("sal")) varexist = true;
  if(ppar->varname == "ptemp" && HaveVar("temp") && HaveVar("sal")) varexist = true;
  if(ppar->varname == "pdens" && (HaveVar("ptemp") || HaveVar("temp")) && HaveVar("sal")) varexist = true;
  if((ppar->varname == "U" || ppar->varname == "U2") && HaveVar("u") && HaveVar("v")) varexist = true;
  if(!varexist) return {nullptr, "Variable " + ppar->varname + " not exists in this dataset"};
 }
 if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger<size_t>();
 if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"};
 real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer);

 if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"};

 {
  auto dom  = DetGeoDomain(lonb, lone);
  real lon1 = ToGeoDomain(args.at("lonb").ToReal(), dom);
  real lon2 = ToGeoDomain(args.at("lone").ToReal(), dom);
  real lat1 = args.at("latb").ToReal();
  real lat2 = args.at("late").ToReal();

  ppar->yb = static_cast<size_t>(Floor((lat1 - latb) / latstep));
  ppar->ye = static_cast<size_t>(Ceil((lat2 - latb) / latstep));
  if(ppar->ye > ny - 1) ppar->ye = ny - 1;
  if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"};

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

  if(ppar->xb == ppar->xe) return {nullptr, "Lonb must be not equal late"};

  if(depth < 0.0 || depth > depths.back())
   ppar->layer = (depth < 0.0) ? 0 : (depths.size() - 1);
  else
   for(size_t i = 0; i < depths.size() - 1; i++)
   {
    if(depth >= depths[i] && depth <= depths[i + 1])
    {
     ppar->layer = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1);
     break;
    }
   }
 }

 pars.SetParameter("variable", ppar->varname);
 pars.SetParameter("depth", Depth(ppar->layer));
 pars.SetParameter("layer", ppar->layer);
 pars.SetParameter("dataset", Title());
 pars.SetParameter("lonb", Lon(ppar->xb));
 pars.SetParameter("latb", Lat(ppar->yb));
 pars.SetParameter("lone", Lon(ppar->xe));
 pars.SetParameter("late", Lat(ppar->ye));

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

LayeredData::Data LayeredData::Read(const BaseParameters* ip, size_t i) const
{
 if(!isOk()) return Data();
 bool nodepth = false;

 auto p               = dynamic_cast<const struct Parameters*>(ip);
 auto [name, id, tid] = VarNameLoc(p->varname, times[i]);
 if(!name.Exist()) // Conversion read
 {
  // ptemp from temp and sal
  if(p->varname == "ptemp")
  {
   struct Parameters params = *p;
   params.varname           = "temp";
   auto temp                = Read(&params, i);
   params.varname           = "sal";
   auto sal                 = Read(&params, i);
   if(!(temp && sal)) return Data();
   auto out = temp;
   for(size_t ind = 0; ind < out.N(); ind++)
   {
    if(temp.IsFill(ind) || sal.IsFill(ind))
     out.V(ind) = out.Fillval();
    else
     out.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
   }
   return out;
  }
  // temp from ptemp and sal
  if(p->varname == "temp")
  {
   struct Parameters params = *p;
   params.varname           = "ptemp";
   auto temp                = Read(&params, i);
   params.varname           = "sal";
   auto sal                 = Read(&params, i);
   if(!(temp && sal)) return Data();
   auto out = temp;
   for(size_t ind = 0; ind < out.N(); ind++)
   {
    if(temp.IsFill(ind) || sal.IsFill(ind))
     out.V(ind) = out.Fillval();
    else
     out.V(ind) = PTemp2Temp(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
   }
   return out;
  }
  // pdens from temp and sal
  if(p->varname == "pdens")
  {
   struct Parameters params    = *p;
   bool              tempispot = HaveVar("ptemp");
   params.varname              = tempispot ? "ptemp" : "temp";
   auto temp                   = Read(&params, i);
   params.varname              = "sal";
   auto sal                    = Read(&params, i);
   if(!(temp && sal)) return Data();
   auto out = temp;
   for(size_t ind = 0; ind < out.N(); ind++)
   {
    if(temp.IsFill(ind) || sal.IsFill(ind))
     out.V(ind) = out.Fillval();
    else
     out.V(ind) = tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind))
                            : Temp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
   }
   return out;
  }
  // U and U2 from u and v
  if(p->varname == "U" || p->varname == "U2")
  {
   bool              square = p->varname == "U2";
   struct Parameters params = *p;
   params.varname           = "u";
   auto u                   = Read(&params, i);
   params.varname           = "v";
   auto v                   = Read(&params, i);
   if(!(u && v)) return Data();
   auto out = u;
   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
 auto head = nc[id]->Header();
 for(const auto& v: head.Variables())
  if(v.Name() == name)
  {
   if(v.Dimensions().size() == 3) nodepth = true;
   if(v.Type().Id() == NC_SHORT) return ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
   if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(nc[id], name, tid, nodepth, p);
  }

 return Data();
}

template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct LayeredData::Parameters* p) const
{
 real     unitmul = 1.0;
 DataType fill;
 real     offset = 0.0, scale = 1.0;

 {
  auto a_fill     = f->A<DataType>(name, "_FillValue");
  auto a_offset_d = f->A<double>(name, "add_offset");
  auto a_scale_d  = f->A<double>(name, "scale_factor");
  auto a_offset_f = f->A<float>(name, "add_offset");
  auto a_scale_f  = f->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 unit = f->A<MString>(name, "units");
 if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0;

 Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep);

 if(p->xb < p->xe)
 {
  auto var = nodepth ? f->V<DataType>(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
                     : f->V<DataType>(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 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 = nodepth ? f->V<DataType>(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
                      : f->V<DataType>(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
  auto var2 = nodepth ? f->V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
                      : f->V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 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;
}