Remap: check range of grid cell center coordinates.

src/Verifygrid.cc

@@ -733,13 +733,13 @@ print_lonlat(int dig, int gridID, size_t gridsize, const Varray<double> &lon, co
   cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_NAME, name, &length);
   const auto xmm = varrayMinMax(gridsize, lon);
   cdoPrint("%10s : %.*g to %.*g degrees", name, dig, xmm.min, dig, xmm.max);
-  if (xmm.min < -720.0 || xmm.max > 720.0) cdoWarning("Longitudes out of bounds!");
+  if (xmm.min < -720.0 || xmm.max > 720.0) cdoWarning("Grid cell center longitudes out of range!");
 
   length = CDI_MAX_NAME;
   cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_NAME, name, &length);
   const auto ymm = varrayMinMax(gridsize, lat);
   cdoPrint("%10s : %.*g to %.*g degrees", name, dig, ymm.min, dig, ymm.max);
-  if (ymm.min < -90.0 || ymm.max > 90.0) cdoWarning("Latitudes out of bounds!");
+  if (ymm.min < -90.0 || ymm.max > 90.0) cdoWarning("Grid cell center latitudes out of range!");
 }
 
 void *

src/remaplib.cc

@@ -197,36 +197,68 @@ remapgrid_get_lonlat(RemapGrid *grid, size_t cell_add)
 }
 
 void
-check_lon_range(size_t nlons, Varray<double> &lons)
+check_lon_range(const char *txt, size_t nlons, Varray<double> &lons)
 {
   assert(!lons.empty());
 
+  if (txt)
+    {
+      double minval = 1.e36;
+      double maxval = -1.e36;
+#ifdef _OPENMP
+#pragma omp parallel for default(none) shared(nlons, lons)  reduction(min : minval) reduction(max : maxval)
+#endif
+      for (size_t i = 0; i < nlons; ++i)
+        {
+          minval = std::min(minval, lons[i]);
+          maxval = std::max(maxval, lons[i]);
+        }
+      if (minval < -PI2 || maxval > 2 * PI2)
+        cdoWarning("%s grid cell center longitudes out of range (min=%g/max=%g)!", txt, RAD2DEG*minval, RAD2DEG*maxval);
+    }
+
 #ifdef _OPENMP
 #pragma omp parallel for default(none) shared(nlons, lons)
 #endif
-  for (size_t n = 0; n < nlons; ++n)
+  for (size_t i = 0; i < nlons; ++i)
     {
       // remove missing values
-      if (lons[n] < -PI2) lons[n] = 0;
-      if (lons[n] > 2 * PI2) lons[n] = PI2;
+      if (lons[i] < -PI2) lons[i] = 0;
+      if (lons[i] > 2 * PI2) lons[i] = PI2;
 
-      if (lons[n] > PI2) lons[n] -= PI2;
-      if (lons[n] < 0.0) lons[n] += PI2;
+      if (lons[i] > PI2) lons[i] -= PI2;
+      if (lons[i] < 0.0) lons[i] += PI2;
     }
 }
 
 void
-check_lat_range(size_t nlats, Varray<double> &lats)
+check_lat_range(const char *txt, size_t nlats, Varray<double> &lats)
 {
   assert(!lats.empty());
 
+  if (txt)
+    {
+      double minval = 1.e36;
+      double maxval = -1.e36;
+#ifdef _OPENMP
+#pragma omp parallel for default(none) shared(nlats, lats)  reduction(min : minval) reduction(max : maxval)
+#endif
+      for (size_t i = 0; i < nlats; ++i)
+        {
+          minval = std::min(minval, lats[i]);
+          maxval = std::max(maxval, lats[i]);
+        }
+      if (minval < -PIH || maxval > PIH)
+        cdoWarning("%s grid cell center latitudes out of range (min=%g/max=%g)!", txt, RAD2DEG*minval, RAD2DEG*maxval);
+    }
+
 #ifdef _OPENMP
 #pragma omp parallel for default(none) shared(nlats, lats)
 #endif
-  for (size_t n = 0; n < nlats; ++n)
+  for (size_t i = 0; i < nlats; ++i)
     {
-      if (lats[n] > PIH) lats[n] = PIH;
-      if (lats[n] < -PIH) lats[n] = -PIH;
+      if (lats[i] > PIH) lats[i] = PIH;
+      if (lats[i] < -PIH) lats[i] = -PIH;
     }
 }
 
@@ -469,16 +501,15 @@ remapDefineGrid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *tx
   if (lgrid_destroy) gridDestroy(gridID);
 
   // Convert longitudes to 0,2pi interval
-  check_lon_range(grid.size, grid.cell_center_lon);
-  if (grid.num_cell_corners && grid.lneed_cell_corners) check_lon_range(grid.num_cell_corners * grid.size, grid.cell_corner_lon);
-
-  //  Make sure input latitude range is within the machine values for +/- pi/2
-  check_lat_range(grid.size, grid.cell_center_lat);
+  check_lon_range(txt, grid.size, grid.cell_center_lon);
+  if (grid.num_cell_corners && grid.lneed_cell_corners) check_lon_range(nullptr, grid.num_cell_corners * grid.size, grid.cell_corner_lon);
 
-  if (grid.num_cell_corners && grid.lneed_cell_corners) check_lat_range(grid.num_cell_corners * grid.size, grid.cell_corner_lat);
+  // Make sure input latitude range is within the machine values for +/- pi/2
+  check_lat_range(txt, grid.size, grid.cell_center_lat);
+  if (grid.num_cell_corners && grid.lneed_cell_corners) check_lat_range(nullptr, grid.num_cell_corners * grid.size, grid.cell_corner_lat);
 }
 
-/*  Compute bounding boxes for restricting future grid searches */
+// Compute bounding boxes for restricting future grid searches
 static void
 cell_bounding_boxes(RemapGrid &grid, float *cell_bound_box, int remap_grid_basis)
 {