Remapstat: added function calc_maxdist().

src/Remapstat.cc

@@ -183,20 +183,30 @@ void read_ybounds_reg2d(int gridID, Varray<double> &ybounds2d)
 static
 double calc_maxdist(size_t i, size_t nv, double plon, double plat, const Varray<double> &xbounds, const Varray<double> &ybounds)
 {
-  cdoAbort("Internal error: calc_maxdist() not implemented!");
+  double p1[3], p2[3];
+  gcLLtoXYZ(plon, plat, p1);
 
   double maxdist = 0.0;
+  for (size_t k = 0; k < nv; ++k)
+    {
+      const auto lons = &xbounds[nv * i];
+      const auto lats = &ybounds[nv * i];
+      gcLLtoXYZ(lons[k], lats[k], p2);
+      const auto sdist = squareDistance(p1, p2);
+      maxdist = std::max(maxdist, sdist);
+    }
+  maxdist = 1.01 * std::sqrt(maxdist);
+
   return maxdist;
 }
 
 static
 double calc_maxdist_rec2d(size_t i, size_t nlon, double plon, double plat, const Varray<double> &xbounds, const Varray<double> &ybounds)
 {
-  constexpr int nv = 4;
-
   double p1[3], p2[3];
   gcLLtoXYZ(plon, plat, p1);
 
+  constexpr int nv = 4;
   double lons[nv], lats[nv];
   const auto iy = i / nlon;
   const auto ix = i - iy * nlon;
@@ -221,6 +231,15 @@ double calc_maxdist_rec2d(size_t i, size_t nlon, double plon, double plat, const
   return maxdist;
 }
 
+static
+size_t find_points(size_t i, size_t nv, size_t nadds, std::vector<size_t> &adds, const Varray<double> &xvals1,
+                   const Varray<double> &yvals1, const Varray<double> &xbounds, const Varray<double> &ybounds)
+{
+  cdoAbort("Internal error: find_points() not implemented!");
+  size_t nvalues = 0;
+  return nvalues;
+}
+
 static
 size_t find_points_rec2d(size_t i, size_t nlon2, size_t nadds, std::vector<size_t> &adds, const Varray<double> &xvals1,
                          const Varray<double> &yvals1, const Varray<double> &xbounds2d, const Varray<double> &ybounds2d)
@@ -321,8 +340,9 @@ Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
       const auto nadds = gridPointSearchDistanceQnearest(gps, maxdist, lon2, lat2, ndist, adds.data(), dist.data());
       // printf("%zu nadds %zu\n", i+1, nadds);
 
-      size_t nvalues = 0;
-      nvalues = find_points_rec2d(i, nlon2, nadds, adds, xvals1, yvals1, xbounds2d, ybounds2d);
+      const auto nvalues = grid2_is_reg2d
+        ? find_points_rec2d(i, nlon2, nadds, adds, xvals1, yvals1, xbounds2d, ybounds2d)
+        : find_points(i, nv, nadds, adds, xvals1, yvals1, xbounds, ybounds);
 
       if (nvalues)
         {