Skip to content

Commit

Permalink
MIR-673 Subsetting for the pan-Arctic regional reanalysis (CARRA2)
Browse files Browse the repository at this point in the history
  • Loading branch information
pmaciel committed Oct 30, 2024
1 parent 49dc995 commit 04cd6c6
Show file tree
Hide file tree
Showing 6 changed files with 90 additions and 19 deletions.
4 changes: 4 additions & 0 deletions etc/mir/grib-input.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
11 changes: 4 additions & 7 deletions src/mir/repres/regular/PolarStereographic.cc
Original file line number Diff line number Diff line change
Expand Up @@ -84,15 +84,12 @@ void PolarStereographic::fillJob(api::MIRJob& job) const {
<< "Nj=" << y()->size() << ";"
<< "grid=" << std::abs(x().step()) << "/" << std::abs(y().step()) << ";"
<< "latitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LAT] << ";"
<< "longitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LON] << ";"
<< "iScansNegatively=" << (x().front() < x().back() ? 0 : 1) << ";"
<< "jScansPositively=" << (y().front() < y().back() ? 1 : 0);
// << "shapeOfTheEarth=" << << ";"
// << "radius=" << << ";"
// << "earthMajorAxis=" << << ";"
// << "earthMinorAxis=" << << ";"
<< "longitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LON];

job.set("grid", grid.str());

// some extra keys are edition-specific, so parent call is here
RegularGrid::fillJob(job);
}


Expand Down
18 changes: 17 additions & 1 deletion src/mir/repres/regular/RegularGrid.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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;
}
Expand Down
2 changes: 2 additions & 0 deletions src/mir/repres/regular/RegularGrid.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
53 changes: 47 additions & 6 deletions src/mir/util/Shape.cc
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,9 @@

#include "mir/util/Shape.h"

#include <sstream>

#include "mir/api/MIRJob.h"
#include "mir/param/MIRParametrisation.h"


Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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
21 changes: 16 additions & 5 deletions src/mir/util/Shape.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand All @@ -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;
Expand Down

0 comments on commit 04cd6c6

Please sign in to comment.