ECC-1818: GRIB Geoiterator issues for Lambert azimuthal equal area

ecmwf · May 2, 2024 · c04bc4f · c04bc4f
1 parent e06656c
commit c04bc4f
Showing 1 changed file with 11 additions and 2 deletions.
diff --git a/src/grib_iterator_class_lambert_azimuthal_equal_area.cc b/src/grib_iterator_class_lambert_azimuthal_equal_area.cc
@@ -209,7 +209,11 @@ static int init_oblate(grib_handle* h,
     sinphi_  = sin(standardParallelInRadians); /*  (P->phi0); */
     Q__sinb1 = pj_qsfn(sinphi_, e, one_es) / Q__qp;
     Q__cosb1 = sqrt(1.0 - Q__sinb1 * Q__sinb1);
-    Q__dd    = cos(standardParallelInRadians) / (sqrt(1. - es * sinphi_ * sinphi_) * Q__rq * Q__cosb1);
+    if (Q__cosb1 == 0) {
+        Q__dd = 1.0;
+    } else {
+        Q__dd = cos(standardParallelInRadians) / (sqrt(1. - es * sinphi_ * sinphi_) * Q__rq * Q__cosb1);
+    }
     Q__ymf   = (Q__xmf = Q__rq) / Q__dd;
     Q__xmf *= Q__dd;
 
@@ -253,7 +257,12 @@ static int init_oblate(grib_handle* h,
                 xy_y *= Q__dd;
                 rho = hypot(xy_x, xy_y);
                 Assert(rho >= EPS10); /* TODO(masn): check */
-                sCe = 2. * asin(.5 * rho / Q__rq);
+                const double asin_arg = (0.5 * rho / Q__rq);
+                if (asin_arg < -1.0 || asin_arg > 1.0) {
+                    grib_context_log(h->context, GRIB_LOG_ERROR, "Invalid value: arcsin argument=%g", asin_arg);
+                    return GRIB_GEOCALCULUS_PROBLEM;
+                }
+                sCe = 2. * asin(asin_arg);
                 cCe = cos(sCe);
                 sCe = sin(sCe);
                 xy_x *= sCe;