diff --git a/etc/mir/grib-input.yaml b/etc/mir/grib-input.yaml index 1d0302bb7..fe6a7bb31 100644 --- a/etc/mir/grib-input.yaml +++ b/etc/mir/grib-input.yaml @@ -19,6 +19,10 @@ md5GridSection=f5dc74ec36353f4c83f7de3bf46e1aef, centre:i=98, gridName=O640, edi md5GridSection=6ae0076681a2541abe053ab52fa60f30, centre:i=84/85, gridType=lambert, edition=1: - first_point_bottom_left: true +# gridType=polar_stereographic does not support the default area method +gridType=polar_stereographic: + - area-mode: mask + # gridType=space_view does not support the default interpolation method gridType=space_view: - interpolation: k-nearest diff --git a/src/mir/repres/regular/PolarStereographic.cc b/src/mir/repres/regular/PolarStereographic.cc index 167d99e00..670f6091b 100644 --- a/src/mir/repres/regular/PolarStereographic.cc +++ b/src/mir/repres/regular/PolarStereographic.cc @@ -12,7 +12,14 @@ #include "mir/repres/regular/PolarStereographic.h" +#include +#include + +#include "mir/api/MIRJob.h" +#include "mir/param/MIRParametrisation.h" +#include "mir/util/Angles.h" #include "mir/util/Exceptions.h" +#include "mir/util/Grib.h" namespace mir::repres::regular { @@ -22,11 +29,67 @@ static const RepresentationBuilder __builder("polar_stereogr PolarStereographic::PolarStereographic(const param::MIRParametrisation& param) : - RegularGrid(param, make_proj_spec(param)) {} + RegularGrid(param, make_proj_spec(param)), + LaDInDegrees_(0), + orientationOfTheGridInDegrees_(0), + southPoleOnProjectionPlane_(false) { + ASSERT(param.get("proj", proj_) && !proj_.empty()); + + param.get("LaDInDegrees", LaDInDegrees_); + param.get("orientationOfTheGridInDegrees", orientationOfTheGridInDegrees_); + + long south = 0; + param.get("southPoleOnProjectionPlane", south); + southPoleOnProjectionPlane_ = south != 0; +} + + +void PolarStereographic::fillGrib(grib_info& info) const { + info.grid.grid_type = CODES_UTIL_GRID_SPEC_POLAR_STEREOGRAPHIC; + info.extra_set("edition", 2L); // write LaDInDegrees + + Point2 first = {x().front(), y().front()}; + Point2 firstLL = grid().projection().lonlat(first); + + info.grid.latitudeOfFirstGridPointInDegrees = firstLL[LLCOORDS::LAT]; + info.grid.longitudeOfFirstGridPointInDegrees = util::normalise_longitude(firstLL[LLCOORDS::LON], 0); + + info.grid.Ni = static_cast(x().size()); + info.grid.Nj = static_cast(y().size()); + + info.grid.iScansNegatively = x().front() < x().back() ? 0L : 1L; + info.grid.jScansPositively = y().front() < y().back() ? 1L : 0L; + + info.extra_set("DxInMetres", std::abs(x().step())); + info.extra_set("DyInMetres", std::abs(y().step())); + info.extra_set("LaDInDegrees", LaDInDegrees_); + info.extra_set("orientationOfTheGridInDegrees", util::normalise_longitude(orientationOfTheGridInDegrees_, 0)); + + // some extra keys are edition-specific, so parent call is here + RegularGrid::fillGrib(info); +} + + +void PolarStereographic::fillJob(api::MIRJob& job) const { + Point2 first = {x().front(), y().front()}; + Point2 firstLL = grid().projection().lonlat(first); + + std::ostringstream grid; + grid << "gridType=polar_stereographic;" + << "proj=" << proj_ << ";" + << "LaDInDegrees=" << LaDInDegrees_ << ";" + << "orientationOfTheGridInDegrees=" << orientationOfTheGridInDegrees_ << ";" + << "southPoleOnProjectionPlane=" << (southPoleOnProjectionPlane_ ? "1" : "0") << ";" + << "Ni=" << x()->size() << ";" + << "Nj=" << y()->size() << ";" + << "grid=" << std::abs(x().step()) << "/" << std::abs(y().step()) << ";" + << "latitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LAT] << ";" + << "longitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LON]; + job.set("grid", grid.str()); -void PolarStereographic::fillGrib(grib_info& /*info*/) const { - NOTIMP; + // some extra keys are edition-specific, so parent call is here + RegularGrid::fillJob(job); } diff --git a/src/mir/repres/regular/PolarStereographic.h b/src/mir/repres/regular/PolarStereographic.h index 35afa7bea..d2f3d8fe7 100644 --- a/src/mir/repres/regular/PolarStereographic.h +++ b/src/mir/repres/regular/PolarStereographic.h @@ -12,6 +12,8 @@ #pragma once +#include + #include "mir/repres/regular/RegularGrid.h" @@ -25,7 +27,8 @@ class PolarStereographic final : public RegularGrid { // -- Constructors - explicit PolarStereographic(const param::MIRParametrisation&); + PolarStereographic(const param::MIRParametrisation&); + PolarStereographic(const PolarStereographic&) = delete; // -- Destructor // None @@ -34,7 +37,8 @@ class PolarStereographic final : public RegularGrid { // None // -- Operators - // None + + PolarStereographic& operator=(const PolarStereographic&) = delete; // -- Methods // None @@ -50,7 +54,11 @@ class PolarStereographic final : public RegularGrid { private: // -- Members - // None + + std::string proj_; + double LaDInDegrees_; + double orientationOfTheGridInDegrees_; + bool southPoleOnProjectionPlane_; // -- Methods // None @@ -58,6 +66,7 @@ class PolarStereographic final : public RegularGrid { // -- Overridden methods void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; // -- Class members // None diff --git a/src/mir/repres/regular/RegularGrid.cc b/src/mir/repres/regular/RegularGrid.cc index e27f10c5f..4f188a0a8 100644 --- a/src/mir/repres/regular/RegularGrid.cc +++ b/src/mir/repres/regular/RegularGrid.cc @@ -23,6 +23,7 @@ #include "eckit/utils/MD5.h" #include "eckit/utils/StringTools.h" +#include "mir/api/MIRJob.h" #include "mir/param/MIRParametrisation.h" #include "mir/repres/Iterator.h" #include "mir/util/Domain.h" @@ -217,11 +218,26 @@ void RegularGrid::fillGrib(grib_info& info) const { shape_.fillGrib(info, grid_.projection().spec()); // scanningMode - info.grid.iScansNegatively = x_.back() < x_.front() ? 1L : 0L; + info.grid.iScansNegatively = x_.front() < x_.back() ? 0L : 1L; info.grid.jScansPositively = y_.front() < y_.back() ? 1L : 0L; } +void RegularGrid::fillJob(api::MIRJob& job) const { + // shape of the reference system + shape_.fillJob(job, grid_.projection().spec()); + + // scanningMode + std::string grid; + ASSERT(job.get("grid", grid) && !grid.empty()); + + grid += ";iScansNegatively=" + std::to_string(x().front() < x().back() ? 0 : 1); + grid += ";jScansPositively=" + std::to_string(y().front() < y().back() ? 1 : 0); + + job.set("grid", grid); +} + + bool RegularGrid::includesNorthPole() const { return bbox_.north() == Latitude::NORTH_POLE; } diff --git a/src/mir/repres/regular/RegularGrid.h b/src/mir/repres/regular/RegularGrid.h index c8d1000d8..ddde26de0 100644 --- a/src/mir/repres/regular/RegularGrid.h +++ b/src/mir/repres/regular/RegularGrid.h @@ -102,7 +102,9 @@ class RegularGrid : public Gridded { bool isPeriodicWestEast() const override; void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; void fillMeshGen(util::MeshGeneratorParameters&) const override; + void reorder(long scanningMode, MIRValuesVector&) const override; void validate(const MIRValuesVector&) const override; void makeName(std::ostream&) const override; diff --git a/src/mir/util/Shape.cc b/src/mir/util/Shape.cc index ea16bee55..f7898ee3b 100644 --- a/src/mir/util/Shape.cc +++ b/src/mir/util/Shape.cc @@ -12,6 +12,9 @@ #include "mir/util/Shape.h" +#include + +#include "mir/api/MIRJob.h" #include "mir/param/MIRParametrisation.h" @@ -50,12 +53,6 @@ Shape::Shape(const Projection::Spec& spec) : edition(0) { } -Shape::Shape(const Shape&) = default; - - -Shape& Shape::operator=(const Shape&) = default; - - void Shape::fillGrib(grib_info& info, const Projection::Spec& spec) const { // GRIB2 encoding of user-provided shape if (edition != 2) { @@ -93,4 +90,48 @@ void Shape::fillGrib(grib_info& info, const Projection::Spec& spec) const { } +void Shape::fillJob(api::MIRJob& job, const Projection::Spec& spec) const { + // GRIB2 encoding of user-provided shape + if (edition != 2) { + return; + } + + std::ostringstream shape; + + // shape given by radius or semi-major/minor axis + if (provided) { + shape << ";shapeOfTheEarth=" << code; + switch (code) { + case 1: + shape << ";radius", spec.getDouble("radius", a); + break; + case 3: + shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis", a) / 1000.; + shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis", b) / 1000.; + break; + case 7: + shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis", a); + shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis", b); + break; + default: + break; + } + } + else if (spec.has("radius")) { + shape << ";shapeOfTheEarth=" << 1L; + shape << ";radius", spec.getDouble("radius"); + } + else if (spec.has("semi_major_axis") && spec.has("semi_minor_axis")) { + shape << ";shapeOfTheEarth=" << 7L; + shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis"); + shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis"); + } + + std::string grid; + ASSERT(job.get("grid", grid) && !grid.empty()); + + job.set("grid", grid + ";" + shape.str()); +} + + } // namespace mir::util diff --git a/src/mir/util/Shape.h b/src/mir/util/Shape.h index a538bd1d6..48937b616 100644 --- a/src/mir/util/Shape.h +++ b/src/mir/util/Shape.h @@ -16,9 +16,14 @@ #include "mir/util/Grib.h" -namespace mir::param { +namespace mir { +namespace param { class MIRParametrisation; -} // namespace mir::param +} +namespace api { +class MIRJob; +} +} // namespace mir namespace mir::util { @@ -30,11 +35,17 @@ struct Shape { Shape(const param::MIRParametrisation&); Shape(const Projection::Spec&); Shape() : code(6L), a(0.), b(0.), provided(false) {} - Shape(const Shape& other); + + Shape(const Shape& other) = default; + Shape(Shape&&) = default; + virtual ~Shape() = default; - Shape& operator=(const Shape& other); - void fillGrib(grib_info& info, const Projection::Spec&) const; + Shape& operator=(const Shape&) = default; + Shape& operator=(Shape&&) = default; + + void fillGrib(grib_info&, const Projection::Spec&) const; + void fillJob(api::MIRJob&, const Projection::Spec&) const; long code; double a; diff --git a/tests/assertions/MIR-673.001.test b/tests/assertions/MIR-673.001.test new file mode 100644 index 000000000..9298d539e --- /dev/null +++ b/tests/assertions/MIR-673.001.test @@ -0,0 +1,10 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--grid=gridType=polar_stereographic +# grib_get +gridType=polar_stereographic +Ni=2869 +Nj=2869 +Dx=2500000 +Dy=2500000 diff --git a/tests/assertions/MIR-673.002.test b/tests/assertions/MIR-673.002.test new file mode 100644 index 000000000..5a2cf8bb1 --- /dev/null +++ b/tests/assertions/MIR-673.002.test @@ -0,0 +1,10 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--area=66.5/-24.5/63.4/-13.5 +# grib_get +gridType=polar_stereographic +Ni=2869 +Nj=2869 +Dx=2500000 +Dy=2500000 diff --git a/tests/assertions/MIR-673.003.fail.test b/tests/assertions/MIR-673.003.fail.test new file mode 100644 index 000000000..d7aef75c2 --- /dev/null +++ b/tests/assertions/MIR-673.003.fail.test @@ -0,0 +1,4 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--area=66.5/-24.5/63.4/-13.5 --area-mode=crop diff --git a/tests/assertions/gridType=polar_stereographic.grib2 b/tests/assertions/gridType=polar_stereographic.grib2 new file mode 120000 index 000000000..6a4ca1f01 --- /dev/null +++ b/tests/assertions/gridType=polar_stereographic.grib2 @@ -0,0 +1 @@ +../data/gridType=polar_stereographic.grib2 \ No newline at end of file diff --git a/tests/data/gridType=polar_stereographic.grib2 b/tests/data/gridType=polar_stereographic.grib2 new file mode 100644 index 000000000..73c9f77e3 Binary files /dev/null and b/tests/data/gridType=polar_stereographic.grib2 differ diff --git a/tests/plans/MIR-673.001.test b/tests/plans/MIR-673.001.test new file mode 100644 index 000000000..abdf907a1 --- /dev/null +++ b/tests/plans/MIR-673.001.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--grid=gridType=polar_stereographic +# plan +Copy[] diff --git a/tests/plans/MIR-673.002.test b/tests/plans/MIR-673.002.test new file mode 100644 index 000000000..dce943e52 --- /dev/null +++ b/tests/plans/MIR-673.002.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--area=66.5/-24.5/63.4/-13.5 +# plan +AreaMasker[bbox=BoundingBox[north=66.5,west=-24.5,south=63.4,east=-13.5]]|Save[output=...] diff --git a/tests/plans/MIR-673.003.test b/tests/plans/MIR-673.003.test new file mode 100644 index 000000000..a8f82df91 --- /dev/null +++ b/tests/plans/MIR-673.003.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--area=66.5/-24.5/63.4/-13.5 --area-mode=crop +# plan +AreaCropper[bbox=BoundingBox[north=66.5,west=-24.5,south=63.4,east=-13.5]]|Save[output=...] diff --git a/tests/plans/MIR-673.004.test b/tests/plans/MIR-673.004.test new file mode 100644 index 000000000..91e5762a0 --- /dev/null +++ b/tests/plans/MIR-673.004.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--same=data.in.MIR-673.003.test +# plan +Gridded2TypedGrid[grid={"LaDInDegrees":"90","Ni":"2869","Nj":"2869","grid":"2500/2500","gridType":"polar_stereographic","iScansNegatively":"0","jScansPositively":"1","latitudeOfFirstGridPointInDegrees":"42.1138","longitudeOfFirstGridPointInDegrees":"-71.6168","orientationOfTheGridInDegrees":"-30","proj":"+proj=stere +lat_ts=90.000003 +lat_0=90 +lon_0=-30.000001 +k_0=1 +x_0=0 +y_0=0 +R=6371229.000000","shapeOfTheEarth":"6","southPoleOnProjectionPlane":"0"},interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/gridType=polar_stereographic.grib2 b/tests/plans/gridType=polar_stereographic.grib2 new file mode 120000 index 000000000..6a4ca1f01 --- /dev/null +++ b/tests/plans/gridType=polar_stereographic.grib2 @@ -0,0 +1 @@ +../data/gridType=polar_stereographic.grib2 \ No newline at end of file