You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
363 lines
12 KiB
363 lines
12 KiB
#pragma once |
|
#include "BFileW.h" |
|
#include "CF.h" |
|
#include "actiondep.h" |
|
#include <memory> |
|
|
|
using michlib::BFileW; |
|
using michlib::Cos; |
|
using michlib::M_PI; |
|
|
|
ADD_ACTION(GenIntFile, genintfile, AdapterSupported<Source> || ReadIsGrid<Source>); |
|
|
|
template<class D> MString ActionGenIntFile_DoAction(const CLArgs& args, D& data); |
|
|
|
template<class D> MString ActionGenIntFile::DoAction(const CLArgs& args, D& data) { return ActionGenIntFile_DoAction(args, data); } |
|
|
|
template<class D> |
|
requires AdapterSupported<D> |
|
MString ActionGenIntFile_DoAction(const CLArgs& args, D& ds) |
|
{ |
|
michlib_internal::ParameterListEx pars; |
|
pars.UsePrefix(""); |
|
pars.SetParameter("source", args.at("source")); |
|
pars.SetParameter("history", args.at("_cmdline")); |
|
|
|
Adapter ad; |
|
{ |
|
auto ret = ds.GetAdapter(args, pars); |
|
if(!ret) return "Error opening source"; |
|
ad = std::move(ret.Value()); |
|
} |
|
|
|
// Modules |
|
{ |
|
auto ret = CF::TEOS10Module(args, pars, ad); |
|
if(!ret) return "Error adding TEOS-10 derived variables"; |
|
} |
|
{ |
|
auto ret = CF::GeostrophicModule(args, pars, ad); |
|
if(!ret) return "Error adding geostrophic velocity"; |
|
} |
|
|
|
if(!args.contains("out")) return "Output file name not specified"; |
|
|
|
MString uvar, vvar; |
|
|
|
// Get U and V |
|
{ |
|
MString upar = args.contains("u") ? args.at("u") : "auto"; |
|
MString vpar = args.contains("v") ? args.at("v") : "auto"; |
|
|
|
if(upar == "auto") |
|
{ |
|
if(ad.HasVar("u")) |
|
uvar = "u"; |
|
else if(ad.HasVar("ugs")) |
|
uvar = "ugs"; |
|
else |
|
return "Can't find u variable"; |
|
} |
|
else |
|
uvar = upar; |
|
|
|
if(vpar == "auto") |
|
{ |
|
MString temp; |
|
for(size_t i = 0; i < uvar.Len(); i++) temp += uvar[i] == 'u' ? "v" : uvar.SubStr(i + 1, 1); |
|
if(temp != uvar && ad.HasVar(temp)) |
|
vvar = temp; |
|
else |
|
return "Can't find v variable"; |
|
} |
|
else |
|
vvar = vpar; |
|
} |
|
|
|
// Check grids |
|
{ |
|
const auto& proj = ad.Vars().at(uvar).proj; |
|
if(proj != ad.Vars().at(vvar).proj) return "U and V variables has different grids"; |
|
if(proj->Type() != Projection::Type::EQC) return "U and V variables has unsupported grid"; |
|
if((ad.Is2D(uvar) && ad.Is3D(vvar)) || (ad.Is3D(uvar) && ad.Is2D(vvar))) return "U and V variables has different type"; |
|
if(ad.Is3D(uvar) && (ad.Vars().at(uvar).vert != ad.Vars().at(vvar).vert)) return "U and V variables has different vertical grids"; |
|
if(ad.Is3D(uvar) && ad.Vars().at(uvar).vert->Nz() != 1) return "U and V variables has multilayer vertical grid"; |
|
} |
|
|
|
// Units |
|
const UtUnit inu(ad.Vars().at(uvar).info->targetunit), inv(ad.Vars().at(vvar).info->targetunit); |
|
const UtUnit cms("cm/s"); |
|
const auto cnvu = Converter(inu, cms); |
|
const auto cnvv = Converter(inv, cms); |
|
|
|
if(!cnvu) return "Can't convert units from \"" + ad.Vars().at(uvar).info->targetunit + "\" to \"cm\""; |
|
if(!cnvv) return "Can't convert units from \"" + ad.Vars().at(vvar).info->targetunit + "\" to \"cm\""; |
|
|
|
const bool needsconvertu = inu != cms; |
|
const bool needsconvertv = inv != cms; |
|
|
|
auto OutU = [needsconvert = needsconvertu, &cnv = cnvu](real r) -> real { return needsconvert ? cnv(r) : r; }; |
|
auto OutV = [needsconvert = needsconvertv, &cnv = cnvv](real r) -> real { return needsconvert ? cnv(r) : r; }; |
|
|
|
const auto& ind = ad.TimeIndexes(); |
|
const auto times = ind | std::views::transform([× = ad.Times()](size_t i) { return times[i]; }); |
|
if(!CF::CheckUniform(times)) return "Time is not uniform"; |
|
|
|
michlib::message("U is \"", uvar, "\", V is \"", vvar, "\""); |
|
pars.SetParameter("U", uvar); |
|
pars.SetParameter("V", vvar); |
|
|
|
const real fillval = 1e10; |
|
|
|
MString name = args.at("out"); |
|
BFileW fw; |
|
fw.Create(name, 2); |
|
fw.SetColumnName(1, "u"); |
|
fw.SetColumnName(2, "v"); |
|
|
|
const auto& proj = ad.Vars().at(uvar).proj; |
|
const size_t Nx = proj->Nx(); |
|
const size_t Ny = proj->Ny(); |
|
{ |
|
real lonb = proj->Lon(0), latb = proj->Lat(0), lone = proj->Lon(Nx - 1), late = proj->Lat(Ny - 1); |
|
fw.UsePrefix(""); |
|
fw.SetParameter("nx", Nx); |
|
fw.SetParameter("ny", Ny); |
|
fw.SetParameter("dx", (lone - lonb) / (Nx - 1) * 60.0); |
|
fw.SetParameter("dy", (late - latb) / (Ny - 1) * 60.0); |
|
fw.SetParameter("nt", times.size()); |
|
fw.SetParameter("dt", (times[1] - times[0]).D()); |
|
fw.SetParameter("FillValue", fillval); |
|
fw.UsePrefix("Info"); |
|
|
|
fw.SetParameter("lonb", lonb); |
|
fw.SetParameter("latb", latb); |
|
fw.SetParameter("lone", lone); |
|
fw.SetParameter("late", late); |
|
fw.SetParameter("DataSource", ad.Title()); |
|
fw.SetParameter("xy2lon", MString(lonb) + "+%x/60"); |
|
fw.SetParameter("xy2lon_rp", MString("POP 60 DIV ") + lonb + " ADD"); |
|
fw.SetParameter("lonlat2x", MString("(%lon-") + lonb + ")*60"); |
|
fw.SetParameter("xy2lat", MString(latb) + "+%y/60"); |
|
fw.SetParameter("xy2lat_rp", MString("EXCH POP 60 DIV ") + latb + " ADD"); |
|
fw.SetParameter("lonlat2y", MString("(%lat-") + latb + ")*60"); |
|
fw.SetParameter("lonlatuv2dotx", "%u*(0.864/1.852)/Cos(%lat*M_PI/180)"); |
|
fw.SetParameter("lonlatuv2doty", "%v*(0.864/1.852)"); |
|
fw.SetParameter("dotxdotylonlat2u_rp", "EXCH POP EXCH POP COSD 1.852 0.864 DIV MUL MUL"); |
|
fw.SetParameter("dotxdotylonlat2v_rp", "POP POP 1.852 0.864 DIV MUL EXCH POP"); |
|
fw.SetParameter("dotxdotylonlat2u", "(1.852/0.864)*Cos(%lat*M_PI/180)*%dotx"); |
|
fw.SetParameter("dotxdotylonlat2v", "(1.852/0.864)*%doty"); |
|
fw.SetParameter("BeginDate", times.front().ToTString()); |
|
fw.SetParameter("EndDate", times.back().ToTString()); |
|
fw.SetParameter("Timestep", (times[1] - times[0]).Seconds()); |
|
fw.SetParameter("DataID", michlib::UniqueId()); |
|
fw.SetParameter("Creation command", args.at("_cmdline")); |
|
fw.SetParameter("U variable", uvar); |
|
fw.SetParameter("V variable", vvar); |
|
} |
|
|
|
for(size_t it = 0; it < times.size(); it++) |
|
{ |
|
michlib::message(times[it].ToTString()); |
|
|
|
const bool is2D = ad.Is2D(uvar); |
|
|
|
std::shared_ptr<Data2D> udata2, vdata2; |
|
std::shared_ptr<Data3D> udata3, vdata3; |
|
|
|
auto U = [&data2 = std::as_const(udata2), &data3 = std::as_const(udata3), is2D = is2D](size_t ix, size_t iy) { return is2D ? (*data2)(ix, iy) : (*data3)(ix, iy, 0); }; |
|
auto V = [&data2 = std::as_const(vdata2), &data3 = std::as_const(vdata3), is2D = is2D](size_t ix, size_t iy) { return is2D ? (*data2)(ix, iy) : (*data3)(ix, iy, 0); }; |
|
|
|
auto Ufill = [&data2 = std::as_const(udata2), &data3 = std::as_const(udata3), is2D = is2D](size_t ix, size_t iy) |
|
{ return is2D ? data2->IsFill(ix, iy) : data3->IsFill(ix, iy, 0); }; |
|
auto Vfill = [&data2 = std::as_const(vdata2), &data3 = std::as_const(vdata3), is2D = is2D](size_t ix, size_t iy) |
|
{ return is2D ? data2->IsFill(ix, iy) : data3->IsFill(ix, iy, 0); }; |
|
|
|
auto IsCoast = [&U = Ufill, &V = Vfill](size_t ix, size_t iy) { return U(ix, iy) || V(ix, iy); }; |
|
|
|
if(is2D) |
|
{ |
|
auto retu = ad.Read2D(uvar, ind[it]); |
|
if(!retu) return "Can't read U data"; |
|
auto retv = ad.Read2D(vvar, ind[it]); |
|
if(!retv) return "Can't read V data"; |
|
udata2 = retu.Value(); |
|
vdata2 = retv.Value(); |
|
} |
|
else |
|
{ |
|
auto retu = ad.Read3D(uvar, ind[it]); |
|
if(!retu) return "Can't read U data"; |
|
auto retv = ad.Read3D(vvar, ind[it]); |
|
if(!retv) return "Can't read V data"; |
|
udata3 = retu.Value(); |
|
vdata3 = retv.Value(); |
|
} |
|
|
|
for(size_t ilat = 0; ilat < Ny; ilat++) |
|
{ |
|
for(size_t ilon = 0; ilon < Nx; ilon++) |
|
{ |
|
if(IsCoast(ilon, ilat)) |
|
{ |
|
fw.Write(fillval); |
|
fw.Write(fillval); |
|
} |
|
else |
|
{ |
|
fw.Write(OutU(U(ilon, ilat)) * (0.864 / 1.852) / Cos(proj->Lat(ilon, ilat) * M_PI / 180.0)); |
|
fw.Write(OutV(V(ilon, ilat)) * (0.864 / 1.852)); |
|
} |
|
} |
|
} |
|
} |
|
|
|
return ""; |
|
} |
|
|
|
template<class D> |
|
requires(!AdapterSupported<D> && ReadIsGrid<D>) |
|
MString ActionGenIntFile_DoAction(const CLArgs& args, D& ds) |
|
{ |
|
auto [reg, regerr] = GetRegion<D>(args); |
|
if(regerr.Exist()) return regerr; |
|
|
|
auto resop = ds.Open(args); |
|
if(resop.Exist()) return "Can't open source: " + resop; |
|
|
|
michlib_internal::ParameterListEx pars; |
|
pars.UsePrefix(""); |
|
pars.SetParameter("source", args.at("source")); |
|
|
|
if(!args.contains("out")) return "Output file name not specified"; |
|
|
|
VelSource mode = VelSource::AUTOSEL; |
|
if(args.contains("geostrophic")) |
|
{ |
|
auto modestr = args.at("geostrophic"); |
|
if(modestr == "" || modestr == "ssh") |
|
mode = VelSource::GEOSTROPHICSSH; |
|
else if(modestr == "surf" || modestr == "surface") |
|
mode = VelSource::GEOSTROPHIC; |
|
else if(modestr == "nongeo") |
|
mode = VelSource::STANDART; |
|
else |
|
return "Parameter \"geostrophic\" have incorrect value " + modestr; |
|
} |
|
|
|
auto [tindexes, err] = GetTIndexes(ds, args, pars); |
|
if(err.Exist()) return err; |
|
if(tindexes.size() < 10) return "Too few time moments"; |
|
// Check isotropy |
|
for(size_t it = 0; it < tindexes.size() - 1; it++) |
|
if(ds.Time(tindexes[it + 1]) - ds.Time(tindexes[it]) != ds.Time(tindexes[1]) - ds.Time(tindexes[0])) return "The time step is not isotropic"; |
|
|
|
std::unique_ptr<const BaseParameters> sourcepars; |
|
if constexpr(ParametersSupported<D>) |
|
{ |
|
if constexpr(ParametersRequiredRegion<D>) |
|
{ |
|
auto [p, err] = ds.Parameters(pars, args, reg); |
|
if(err.Exist()) return err; |
|
sourcepars.reset(p); |
|
} |
|
else |
|
{ |
|
auto [p, err] = ds.Parameters(pars, args); |
|
if(err.Exist()) return err; |
|
sourcepars.reset(p); |
|
} |
|
} |
|
auto p = sourcepars.get(); |
|
|
|
MString name = args.at("out"); |
|
BFileW fw; |
|
fw.Create(name, 2); |
|
fw.SetColumnName(1, "u"); |
|
fw.SetColumnName(2, "v"); |
|
|
|
for(size_t it = 0; it < tindexes.size(); it++) |
|
{ |
|
auto data = ReadUV(ds, p, tindexes[it], mode); |
|
if(!data) return "Can't read data"; |
|
if(it == 0) |
|
{ |
|
real lonb = data.Lon(0, 0), latb = data.Lat(0, 0), lone = data.Lon(data.Nx() - 1, data.Ny() - 1), late = data.Lat(data.Nx() - 1, data.Ny() - 1); |
|
fw.UsePrefix(""); |
|
fw.SetParameter("nx", data.Nx()); |
|
fw.SetParameter("ny", data.Ny()); |
|
fw.SetParameter("dx", data.XStep() * 60.0); |
|
fw.SetParameter("dy", data.YStep() * 60.0); |
|
fw.SetParameter("nt", tindexes.size()); |
|
fw.SetParameter("dt", (ds.Time(tindexes[1]) - ds.Time(tindexes[0])).D()); |
|
fw.SetParameter("FillValue", data.Fillval()); |
|
fw.UsePrefix("Info"); |
|
switch(mode) |
|
{ |
|
case(VelSource::AUTOSEL): |
|
{ |
|
fw.SetParameter("Mode", "autoselect between standart and geostrophic"); |
|
break; |
|
} |
|
case(VelSource::STANDART): |
|
{ |
|
fw.SetParameter("Mode", "standart"); |
|
break; |
|
} |
|
case(VelSource::GEOSTROPHIC): |
|
{ |
|
fw.SetParameter("Mode", "geostrophic"); |
|
break; |
|
} |
|
case(VelSource::GEOSTROPHICSSH): |
|
{ |
|
fw.SetParameter("Mode", "geostrophic from ssh"); |
|
break; |
|
} |
|
} |
|
fw.SetParameter("lonb", lonb); |
|
fw.SetParameter("latb", latb); |
|
fw.SetParameter("lone", lone); |
|
fw.SetParameter("late", late); |
|
fw.SetParameter("DataSource", ds.DataTitle()); |
|
fw.SetParameter("xy2lon", MString(lonb) + "+%x/60"); |
|
fw.SetParameter("xy2lon_rp", MString("POP 60 DIV ") + lonb + " ADD"); |
|
fw.SetParameter("lonlat2x", MString("(%lon-") + lonb + ")*60"); |
|
fw.SetParameter("xy2lat", MString(latb) + "+%y/60"); |
|
fw.SetParameter("xy2lat_rp", MString("EXCH POP 60 DIV ") + latb + " ADD"); |
|
fw.SetParameter("lonlat2y", MString("(%lat-") + latb + ")*60"); |
|
fw.SetParameter("lonlatuv2dotx", "%u*(0.864/1.852)/Cos(%lat*M_PI/180)"); |
|
fw.SetParameter("lonlatuv2doty", "%v*(0.864/1.852)"); |
|
fw.SetParameter("dotxdotylonlat2u_rp", "EXCH POP EXCH POP COSD 1.852 0.864 DIV MUL MUL"); |
|
fw.SetParameter("dotxdotylonlat2v_rp", "POP POP 1.852 0.864 DIV MUL EXCH POP"); |
|
fw.SetParameter("dotxdotylonlat2u", "(1.852/0.864)*Cos(%lat*M_PI/180)*%dotx"); |
|
fw.SetParameter("dotxdotylonlat2v", "(1.852/0.864)*%doty"); |
|
fw.SetParameter("BeginDate", ds.Time(tindexes.front()).ToTString()); |
|
fw.SetParameter("EndDate", ds.Time(tindexes.back()).ToTString()); |
|
fw.SetParameter("Timestep", (ds.Time(tindexes[1]) - ds.Time(tindexes[0])).Seconds()); |
|
fw.SetParameter("DataID", michlib::UniqueId()); |
|
fw.SetParameter("Creation command", args.at("_cmdline")); |
|
} |
|
|
|
for(size_t ilat = 0; ilat < data.Ny(); ilat++) |
|
{ |
|
for(size_t ilon = 0; ilon < data.Nx(); ilon++) |
|
{ |
|
if(data.IsCoast(ilon, ilat)) |
|
{ |
|
fw.Write(data.Fillval()); |
|
fw.Write(data.Fillval()); |
|
} |
|
else |
|
{ |
|
fw.Write(data.U(ilon, ilat) * (0.864 / 1.852) / Cos(data.Lat(ilon, ilat) * M_PI / 180.0)); |
|
fw.Write(data.V(ilon, ilat) * (0.864 / 1.852)); |
|
} |
|
} |
|
} |
|
} |
|
fw.Finalize(); |
|
fw.Close(); |
|
|
|
return ""; |
|
};
|
|
|