diff --git a/include/HYCOM.h b/include/HYCOM.h index aec585c..3bc62fd 100644 --- a/include/HYCOM.h +++ b/include/HYCOM.h @@ -11,13 +11,14 @@ #include 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 urls; std::vector depths; std::vector times; + MString lonname, latname; size_t nx, ny; real lonb, latb, lone, late; real lonstep, latstep; @@ -56,7 +58,7 @@ class HYCOMData real offset = 0.0, scale = 1.0; { - auto a_fill = f.A(name, "_FillValue"); + auto a_fill = f.A(name, "_FillValue"); auto a_offset_d = f.A(name, "add_offset"); auto a_scale_d = f.A(name, "scale_factor"); auto a_offset_f = f.A(name, "add_offset"); @@ -76,8 +78,8 @@ class HYCOMData if(p->xb < p->xe) { - auto var = nodepth ? f.V(name, {"lon", p->xb, p->xe - p->xb + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(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(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(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(name, {"lon", p->xb}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(name, {"lon", p->xb}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); - auto var2 = nodepth ? f.V(name, {"lon", 0, p->xe + 1}, {"lat", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(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(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); + auto var2 = nodepth ? f.V(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(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("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(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("lon"); auto lats = nc[0].V("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(); diff --git a/include/NEMO.h b/include/NEMO.h index 949f3b9..447d8ce 100644 --- a/include/NEMO.h +++ b/include/NEMO.h @@ -11,13 +11,14 @@ #include 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 nc; + std::vector urls; std::vector depths; std::vector times; + MString lonname, latname; size_t nx, ny; real lonb, latb, lone, late; real lonstep, latstep; @@ -91,13 +94,17 @@ class NEMOData real offset = 0.0, scale = 1.0; { - auto a_fill = f.A(name, "_FillValue"); - auto a_offset = f.A(name, "add_offset"); - auto a_scale = f.A(name, "scale_factor"); + auto a_fill = f.A(name, "_FillValue"); + auto a_offset_d = f.A(name, "add_offset"); + auto a_scale_d = f.A(name, "scale_factor"); + auto a_offset_f = f.A(name, "add_offset"); + auto a_scale_f = f.A(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(name, "units"); @@ -107,8 +114,8 @@ class NEMOData if(p->xb < p->xe) { - auto var = nodepth ? f.V(name, {"longitude", p->xb, p->xe - p->xb + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(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(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(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(name, {"longitude", p->xb}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(name, {"longitude", p->xb}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); - auto var2 = nodepth ? f.V(name, {"longitude", 0, p->xe + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}) - : f.V(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(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1}); + auto var2 = nodepth ? f.V(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) + : f.V(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 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("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(timeD ? timeD(i) : timeF(i)) * 3600; - auto nlon = nc[0].D("longitude"); - auto nlat = nc[0].D("latitude"); - if(!(nlon && nlat)) + MDateTime refdate; { - nc.clear(); - return "Can't get number of longitudes/latitudes"; + auto units = nc[0].Attribute("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 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(timeD ? timeD(i) : timeF(i)) * 3600; auto lons = nc[0].V("longitude"); auto lats = nc[0].V("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();