Projection base implementation speedup on calculation of bounding box…

…es; MercatorProjection fix on handling extrema
ecmwf · Apr 2, 2024 · 3f2e234 · 3f2e234
1 parent 5c1b201
commit 3f2e234
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 2 deletions.
diff --git a/src/atlas/projection/detail/MercatorProjection.cc b/src/atlas/projection/detail/MercatorProjection.cc
@@ -10,9 +10,11 @@
 
 #include <cmath>
 #include <functional>
+#include <limits>
 #include <sstream>
 
 #include "eckit/config/Parametrisation.h"
+#include "eckit/types/FloatCompare.h"
 #include "eckit/utils/Hash.h"
 
 #include "atlas/projection/detail/MercatorProjection.h"
@@ -96,6 +98,18 @@ void MercatorProjectionT<Rotation>::lonlat2xy(double crd[]) const {
     rotation_.unrotate(crd);
 
     // then project
+    if (eckit::types::is_approximately_equal<double>(crd[LAT], 90., 1e-3)) {
+        crd[XX] = false_easting_;
+        crd[YY] = std::numeric_limits<double>::infinity();
+        return;
+    }
+
+    if (eckit::types::is_approximately_equal<double>(crd[LAT], -90., 1e-3)) {
+        crd[XX] = false_easting_;
+        crd[YY] = -std::numeric_limits<double>::infinity();
+        return;
+    }
+
     crd[XX] = k_radius_ * (D2R(normalise_mercator_(crd[LON]) - lon0_));
     crd[YY] = k_radius_ * 0.5 * std::log(t(crd[LAT]));
     crd[XX] += false_easting_;

diff --git a/src/atlas/projection/detail/ProjectionImpl.cc b/src/atlas/projection/detail/ProjectionImpl.cc
@@ -261,8 +261,8 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain)
     // 2. locate latitude extrema by checking if poles are included (in the un-projected frame) and if not, find extrema
     // not at the corners by refining iteratively
 
-    bounds.includesNorthPole(bounds.includesNorthPole() || rect.contains(xy({0, 90})));
-    bounds.includesSouthPole(bounds.includesSouthPole() || rect.contains(xy({0, -90})));
+    bounds.includesNorthPole(bounds.includesNorthPole() || rect.contains(xy(PointLonLat{0, 90 - h_deg})));
+    bounds.includesSouthPole(bounds.includesSouthPole() || rect.contains(xy(PointLonLat{0, -90 + h_deg})));
 
     for (auto [A, B] : segments) {
         if (!bounds.includesNorthPole() || !bounds.includesSouthPole()) {