Browse Source

Prepare to merge HYCOM and NEMO classes

interpolate
Michael Uleysky 2 years ago
parent
commit
b039db5beb
  1. 54
      include/HYCOM.h
  2. 89
      include/NEMO.h

54
include/HYCOM.h

@ -11,13 +11,14 @@
#include <sys/types.h>
using michlib::Ceil;
using michlib::DetGeoDomain;
using michlib::errmessage;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;
using michlib::MDateTime;
using michlib::NCFileA;
using michlib::ToGeoDomainPos;
using michlib::ToGeoDomain;
class HYCOMData
{
@ -40,6 +41,7 @@ class HYCOMData
std::vector<MString> urls;
std::vector<real> depths;
std::vector<MDateTime> times;
MString lonname, latname;
size_t nx, ny;
real lonb, latb, lone, late;
real lonstep, latstep;
@ -76,8 +78,8 @@ class HYCOMData
if(p->xb < p->xe)
{
auto var = nodepth ? f.V<DataType>(name, {"lon", p->xb, p->xe - p->xb + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"lon", p->xb, p->xe - p->xb + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
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();
@ -91,10 +93,10 @@ class HYCOMData
}
else
{
auto var1 = nodepth ? f.V<DataType>(name, {"lon", p->xb}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"lon", p->xb}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
auto var2 = nodepth ? f.V<DataType>(name, {"lon", 0, p->xe + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"lon", 0, p->xe + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
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++)
@ -189,7 +191,7 @@ class HYCOMData
nc.resize(urls.size());
for(size_t i = 0; i < urls.size(); i++)
{
nc[i].Reset(urls[i]+"#cache&noprefetch");
nc[i].Reset(urls[i] + "#cache&noprefetch");
if(!nc[i])
{
nc.clear();
@ -197,6 +199,27 @@ class HYCOMData
}
}
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";
}
auto rdepths = nc[0].V<double>("depth");
if(!rdepths)
{
@ -240,16 +263,6 @@ class HYCOMData
timeD ? times.resize(timeD.DimLen(0)) : times.resize(timeF.DimLen(0));
for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(timeD ? timeD(i) : timeF(i)) * 3600;
auto nlon = nc[0].D("lon");
auto nlat = nc[0].D("lat");
if(!(nlon && nlat))
{
nc.clear();
return "Can't get number of longitudes/latitudes";
}
nx = nlon;
ny = nlat;
auto lons = nc[0].V<double>("lon");
auto lats = nc[0].V<double>("lat");
if(!(lons && lats))
@ -322,8 +335,9 @@ class HYCOMData
if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"};
{
real lon1 = ToGeoDomainPos(args.at("lonb").ToReal());
real lon2 = ToGeoDomainPos(args.at("lone").ToReal());
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();

89
include/NEMO.h

@ -11,13 +11,14 @@
#include <sys/types.h>
using michlib::Ceil;
using michlib::DetGeoDomain;
using michlib::errmessage;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;
using michlib::MDateTime;
using michlib::NCFileA;
using michlib::ToGeoDomainNeg;
using michlib::ToGeoDomain;
class NEMOData
{
@ -37,8 +38,10 @@ class NEMOData
};
std::vector<NCFileA> nc;
std::vector<MString> urls;
std::vector<real> depths;
std::vector<MDateTime> times;
MString lonname, latname;
size_t nx, ny;
real lonb, latb, lone, late;
real lonstep, latstep;
@ -92,12 +95,16 @@ class NEMOData
{
auto a_fill = f.A<DataType>(name, "_FillValue");
auto a_offset = f.A<double>(name, "add_offset");
auto a_scale = f.A<double>(name, "scale_factor");
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) offset = a_offset;
if(a_scale) scale = a_scale;
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");
@ -107,8 +114,8 @@ class NEMOData
if(p->xb < p->xe)
{
auto var = nodepth ? f.V<DataType>(name, {"longitude", p->xb, p->xe - p->xb + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"longitude", p->xb, p->xe - p->xb + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
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();
@ -121,10 +128,10 @@ class NEMOData
}
else
{
auto var1 = nodepth ? f.V<DataType>(name, {"longitude", p->xb}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"longitude", p->xb}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
auto var2 = nodepth ? f.V<DataType>(name, {"longitude", 0, p->xe + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"longitude", 0, p->xe + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
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++)
@ -205,7 +212,6 @@ class NEMOData
nc.clear();
if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl);
std::vector<MString> urls;
if(dataset == "DT")
{
urls = {"https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m"};
@ -233,7 +239,7 @@ class NEMOData
nc.resize(urls.size());
for(size_t i = 0; i < urls.size(); i++)
{
nc[i].Reset(urls[i]);
nc[i].Reset(urls[i] + "#cache&noprefetch");
if(!nc[i])
{
nc.clear();
@ -241,6 +247,27 @@ class NEMOData
}
}
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";
}
auto rdepths = nc[0].V<float>("depth");
if(!rdepths)
{
@ -257,19 +284,32 @@ class NEMOData
nc.clear();
return "Can't read times";
}
MDateTime refdate("1950-01-01");
timeD ? times.resize(timeD.DimLen(0)) : times.resize(timeF.DimLen(0));
for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(timeD ? timeD(i) : timeF(i)) * 3600;
auto nlon = nc[0].D("longitude");
auto nlat = nc[0].D("latitude");
if(!(nlon && nlat))
MDateTime refdate;
{
auto units = nc[0].Attribute<MString>("time", "units");
if(!units)
{
nc.clear();
return "Can't read refdate";
}
MString rstr;
auto words = michlib::Split_on_words(units);
auto ci = words.begin();
if(ci != words.end()) ci++; // skip "hours"
if(ci != words.end()) ci++; // skip "since"
if(ci != words.end()) rstr = *ci; // Day
if(ci != words.end()) ci++;
if(ci != words.end()) rstr += " " + *ci; // Hours
if(!refdate.FromString(rstr))
{
nc.clear();
return "Can't get number of longitudes/latitudes";
return "Can't parse " + rstr + " to refdate";
}
}
nx = nlon;
ny = nlat;
timeD ? times.resize(timeD.DimLen(0)) : times.resize(timeF.DimLen(0));
for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(timeD ? timeD(i) : timeF(i)) * 3600;
auto lons = nc[0].V<float>("longitude");
auto lats = nc[0].V<float>("latitude");
@ -345,8 +385,9 @@ class NEMOData
if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"};
{
real lon1 = ToGeoDomainNeg(args.at("lonb").ToReal());
real lon2 = ToGeoDomainNeg(args.at("lone").ToReal());
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();

Loading…
Cancel
Save