Added 2nd version of gridPointSearchCreate().

src/grid_point_search.cc

@@ -37,7 +37,7 @@ extern "C"
 constexpr double PI = M_PI;
 constexpr double PI2 = 2.0 * PI;
 
-PointSearchMethod pointSearchMethod(PointSearchMethod::nanoflann);
+PointSearchMethod pointSearchMethod(PointSearchMethod::undefined);
 
 struct gpsFull
 {
@@ -88,11 +88,8 @@ struct PointCloud
   bool
   kdtree_get_bbox(BBOX &bb) const
   {
-    for (unsigned j = 0; j < 3; ++j)
-      {
-        bb[j].low = min[j];
-        bb[j].high = max[j];
-      }
+    for (unsigned i = 0; i < 3; ++i) bb[i].low = min[i];
+    for (unsigned i = 0; i < 3; ++i) bb[i].high = max[i];
     return true;
   }
   // bool kdtree_get_bbox(BBOX& /* bb */) const { return false; }
@@ -102,7 +99,7 @@ using nfTree_t
     = nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, PointCloud<double>>, PointCloud<double>, 3>;
 
 static double
-cdoDefaultDefaultRadius(void)
+cdoDefaultRadius(void)
 {
   extern double pointSearchRadius;
 
@@ -176,23 +173,35 @@ gridPointSearchCreateReg2d(GridPointSearch &gps, bool xIsCyclic, size_t dims[2],
       sinlat[n] = std::sin(lats[n]);
     }
 
-  gps.searchRadius = cdoDefaultDefaultRadius();
+  gps.searchRadius = cdoDefaultRadius();
 
   gps.in_use = true;
 }
 
+static inline void
+min_point(double *min, double *point)
+{
+  for (unsigned i = 0; i < 3; ++i) min[i] = (point[i] < min[i]) ? point[i] : min[i];
+}
+
+static inline void
+max_point(double *max, double *point)
+{
+  for (unsigned i = 0; i < 3; ++i) max[i] = (point[i] > max[i]) ? point[i] : max[i];
+}
+
 template <typename T>
-void
+static void
 adjust_bbox_min(T *min)
 {
-  for (unsigned j = 0; j < 3; ++j) min[j] = min[j] < 0 ? min[j] * 1.001 : min[j] * 0.999;
+  for (unsigned i = 0; i < 3; ++i) min[i] = (min[i] < 0) ? min[i] * 1.001 : min[i] * 0.999;
 }
 
 template <typename T>
-void
+static void
 adjust_bbox_max(T *max)
 {
-  for (unsigned j = 0; j < 3; ++j) max[j] = max[j] < 0 ? max[j] * 0.999 : max[j] * 1.001;
+  for (unsigned i = 0; i < 3; ++i) max[i] = (max[i] < 0) ? max[i] * 0.999 : max[i] * 1.001;
 }
 
 static void *
@@ -211,18 +220,15 @@ gps_create_kdtree(size_t n, const Varray<double> &lons, const Varray<double> &la
     {
       auto &point = pointlist[i].point;
       gcLLtoXYZ(lons[i], lats[i], point);
-      for (unsigned j = 0; j < 3; ++j)
-        {
-          min[j] = point[j] < min[j] ? point[j] : min[j];
-          max[j] = point[j] > max[j] ? point[j] : max[j];
-        }
+      min_point(min, point);
+      max_point(max, point);
       pointlist[i].index = i;
     }
 
   adjust_bbox_min(min);
   adjust_bbox_max(max);
-  for (unsigned j = 0; j < 3; ++j) gps.min[j] = min[j];
-  for (unsigned j = 0; j < 3; ++j) gps.max[j] = max[j];
+  for (unsigned i = 0; i < 3; ++i) gps.min[i] = min[i];
+  for (unsigned i = 0; i < 3; ++i) gps.max[i] = max[i];
 
   if (Options::cdoVerbose) cdoPrint("BBOX: min=%g/%g/%g  max=%g/%g/%g", min[0], min[1], min[2], max[0], max[1], max[2]);
 
@@ -252,22 +258,19 @@ gps_create_nanoflann(size_t n, const Varray<double> &lons, const Varray<double>
       pointcloud->pts[i].x = point[0];
       pointcloud->pts[i].y = point[1];
       pointcloud->pts[i].z = point[2];
-      for (unsigned j = 0; j < 3; ++j)
-        {
-          min[j] = point[j] < min[j] ? point[j] : min[j];
-          max[j] = point[j] > max[j] ? point[j] : max[j];
-        }
+      min_point(min, point);
+      max_point(max, point);
     }
 
   gps.pointcloud = (void *) pointcloud;
 
   adjust_bbox_min(min);
   adjust_bbox_max(max);
-  for (unsigned j = 0; j < 3; ++j) gps.min[j] = min[j];
-  for (unsigned j = 0; j < 3; ++j) gps.max[j] = max[j];
+  for (unsigned i = 0; i < 3; ++i) gps.min[i] = min[i];
+  for (unsigned i = 0; i < 3; ++i) gps.max[i] = max[i];
 
-  for (unsigned j = 0; j < 3; ++j) pointcloud->min[j] = min[j];
-  for (unsigned j = 0; j < 3; ++j) pointcloud->max[j] = max[j];
+  for (unsigned i = 0; i < 3; ++i) pointcloud->min[i] = min[i];
+  for (unsigned i = 0; i < 3; ++i) pointcloud->max[i] = max[i];
 
   // construct a kd-tree index:
   nfTree_t *nft = new nfTree_t(3 /*dim*/, *pointcloud, nanoflann::KDTreeSingleIndexAdaptorParams(50 /* max leaf */));
@@ -293,17 +296,14 @@ gps_create_spherepart(size_t n, const Varray<double> &lons, const Varray<double>
     {
       double *point = coordinates_xyz[i];
       gcLLtoXYZ(lons[i], lats[i], point);
-      for (unsigned j = 0; j < 3; ++j)
-        {
-          min[j] = point[j] < min[j] ? point[j] : min[j];
-          max[j] = point[j] > max[j] ? point[j] : max[j];
-        }
+      min_point(min, point);
+      max_point(max, point);
     }
 
   adjust_bbox_min(min);
   adjust_bbox_max(max);
-  for (unsigned j = 0; j < 3; ++j) gps.min[j] = min[j];
-  for (unsigned j = 0; j < 3; ++j) gps.max[j] = max[j];
+  for (unsigned i = 0; i < 3; ++i) gps.min[i] = min[i];
+  for (unsigned i = 0; i < 3; ++i) gps.max[i] = max[i];
 
   if (Options::cdoVerbose) cdoPrint("BBOX: min=%g/%g/%g  max=%g/%g/%g", min[0], min[1], min[2], max[0], max[1], max[2]);
 
@@ -328,12 +328,8 @@ gps_destroy_full(void *search_container)
   auto full = (gpsFull *) search_container;
   if (full)
     {
-      if (full->pts)
-        {
-          free(full->pts[0]);
-          free(full->pts);
-        }
-
+      if (full->pts) free(full->pts[0]);
+      if (full->pts) free(full->pts);
       free(full);
     }
 }
@@ -371,6 +367,52 @@ gps_create_full(size_t n, const Varray<double> &lons, const Varray<double> &lats
   return (void *) full;
 }
 
+static void
+print_method(PointSearchMethod method)
+{
+  // clang-format off
+  if      (method == PointSearchMethod::kdtree)     cdoPrint("Point search method: kdtree");
+  else if (method == PointSearchMethod::full)       cdoPrint("Point search method: full");
+  else if (method == PointSearchMethod::nanoflann)  cdoPrint("Point search method: nanoflann");
+  else if (method == PointSearchMethod::spherepart) cdoPrint("Point search method: spherepart");
+  // clang-format on
+}
+
+void
+gridPointSearchCreate(GridPointSearch &gps, const Varray<double> &lons, const Varray<double> &lats, PointSearchMethod method)
+{
+  auto n = lons.size();
+
+  gps.is_cyclic = false;
+  gps.is_curve = false;
+  gps.dims[0] = n;
+  gps.dims[1] = 0;
+
+  gps.method = method;
+  if (pointSearchMethod != PointSearchMethod::undefined) gps.method = pointSearchMethod;
+  if (gps.method == PointSearchMethod::latbins) gps.method = PointSearchMethod::spherepart;
+
+  gps.n = n;
+  if (n == 0) return;
+
+  gps.plons = lons.data();
+  gps.plats = lats.data();
+
+  if (Options::cdoVerbose) print_method(gps.method);
+
+  // clang-format off
+  if      (gps.method == PointSearchMethod::kdtree)     gps.search_container = gps_create_kdtree(n, lons, lats, gps);
+  else if (gps.method == PointSearchMethod::full)       gps.search_container = gps_create_full(n, lons, lats);
+  else if (gps.method == PointSearchMethod::nanoflann)  gps.search_container = gps_create_nanoflann(n, lons, lats, gps);
+  else if (gps.method == PointSearchMethod::spherepart) gps.search_container = gps_create_spherepart(n, lons, lats, gps);
+  else cdoAbort("%s::method undefined!", __func__);
+  // clang-format on
+
+  gps.searchRadius = cdoDefaultRadius();
+
+  gps.in_use = true;
+}
+
 void
 gridPointSearchCreate(GridPointSearch &gps, bool xIsCyclic, size_t dims[2], size_t n, const Varray<double> &lons,
                       const Varray<double> &lats)
@@ -380,7 +422,7 @@ gridPointSearchCreate(GridPointSearch &gps, bool xIsCyclic, size_t dims[2], size
   gps.dims[0] = dims[0];
   gps.dims[1] = dims[1];
 
-  gps.method = pointSearchMethod;
+  if (pointSearchMethod != PointSearchMethod::undefined) gps.method = pointSearchMethod;
   if (gps.method == PointSearchMethod::latbins) gps.method = PointSearchMethod::nanoflann;
 
   gps.n = n;
@@ -389,15 +431,9 @@ gridPointSearchCreate(GridPointSearch &gps, bool xIsCyclic, size_t dims[2], size
   gps.plons = lons.data();
   gps.plats = lats.data();
 
-  // clang-format off
-  if (Options::cdoVerbose)
-    {
-      if      (gps.method == PointSearchMethod::kdtree)     cdoPrint("Point search method: kdtree");
-      else if (gps.method == PointSearchMethod::full)       cdoPrint("Point search method: full");
-      else if (gps.method == PointSearchMethod::nanoflann)  cdoPrint("Point search method: nanoflann");
-      else if (gps.method == PointSearchMethod::spherepart) cdoPrint("Point search method: spherepart");
-    }
+  if (Options::cdoVerbose) print_method(gps.method);
 
+  // clang-format off
   if      (gps.method == PointSearchMethod::kdtree)     gps.search_container = gps_create_kdtree(n, lons, lats, gps);
   else if (gps.method == PointSearchMethod::full)       gps.search_container = gps_create_full(n, lons, lats);
   else if (gps.method == PointSearchMethod::nanoflann)  gps.search_container = gps_create_nanoflann(n, lons, lats, gps);
@@ -405,7 +441,7 @@ gridPointSearchCreate(GridPointSearch &gps, bool xIsCyclic, size_t dims[2], size
   else cdoAbort("%s::method undefined!", __func__);
   // clang-format on
 
-  gps.searchRadius = cdoDefaultDefaultRadius();
+  gps.searchRadius = cdoDefaultRadius();
 
   gps.in_use = true;
 }
@@ -457,8 +493,8 @@ gps_nearest_kdtree(void *search_container, double lon, double lat, double search
   gcLLtoXYZ(lon, lat, query_pt);
 
   if (!gps.extrapolate)
-    for (unsigned j = 0; j < 3; ++j)
-      if (query_pt[j] < gps.min[j] || query_pt[j] > gps.max[j]) return 0;
+    for (unsigned i = 0; i < 3; ++i)
+      if (query_pt[i] < gps.min[i] || query_pt[i] > gps.max[i]) return 0;
 
   const auto node = kd_nearest(kdt->node, query_pt, &sqrDist, 3);
 
@@ -485,8 +521,8 @@ gps_nearest_nanoflann(void *search_container, double lon, double lat, double sea
   gcLLtoXYZ(lon, lat, query_pt);
 
   if (!gps.extrapolate)
-    for (unsigned j = 0; j < 3; ++j)
-      if (query_pt[j] < gps.min[j] || query_pt[j] > gps.max[j]) return 0;
+    for (unsigned i = 0; i < 3; ++i)
+      if (query_pt[i] < gps.min[i] || query_pt[i] > gps.max[i]) return 0;
 
   const size_t num_results = 1;
   size_t retIndex;
@@ -513,8 +549,8 @@ gps_nearest_spherepart(void *search_container, double lon, double lat, double se
   gcLLtoXYZ(lon, lat, query_pt[0]);
 
   if (!gps.extrapolate)
-    for (unsigned j = 0; j < 3; ++j)
-      if (query_pt[0][j] < gps.min[j] || query_pt[0][j] > gps.max[j]) return 0;
+    for (unsigned i = 0; i < 3; ++i)
+      if (query_pt[0][i] < gps.min[i] || query_pt[0][i] > gps.max[i]) return 0;
 
   size_t local_point_ids_array_size = 0;
   size_t num_local_point_ids;
@@ -648,8 +684,8 @@ gps_qnearest_kdtree(GridPointSearch &gps, double lon, double lat, double searchR
   gcLLtoXYZ(lon, lat, query_pt);
 
   if (!gps.extrapolate)
-    for (unsigned j = 0; j < 3; ++j)
-      if (query_pt[j] < gps.min[j] || query_pt[j] > gps.max[j]) return nadds;
+    for (unsigned i = 0; i < 3; ++i)
+      if (query_pt[i] < gps.min[i] || query_pt[i] > gps.max[i]) return nadds;
 
   if (gps.in_use)
     {
@@ -694,8 +730,8 @@ gps_qnearest_nanoflann(const GridPointSearch &gps, double lon, double lat, doubl
   gcLLtoXYZ(lon, lat, query_pt);
 
   if (!gps.extrapolate)
-    for (unsigned j = 0; j < 3; ++j)
-      if (query_pt[j] < gps.min[j] || query_pt[j] > gps.max[j]) return nadds;
+    for (unsigned i = 0; i < 3; ++i)
+      if (query_pt[i] < gps.min[i] || query_pt[i] > gps.max[i]) return nadds;
 
   nadds = nft->knnRangeSearch(&query_pt[0], sqrDistMax, nnn, &adds[0], &dist[0]);
   for (size_t i = 0; i < nadds; ++i) dist[i] = std::sqrt(dist[i]);
@@ -714,8 +750,8 @@ gps_qnearest_spherepart(GridPointSearch &gps, double lon, double lat, double sea
       gcLLtoXYZ(lon, lat, query_pt);
 
       if (!gps.extrapolate)
-        for (unsigned j = 0; j < 3; ++j)
-          if (query_pt[j] < gps.min[j] || query_pt[j] > gps.max[j]) return nadds;
+        for (unsigned i = 0; i < 3; ++i)
+          if (query_pt[i] < gps.min[i] || query_pt[i] > gps.max[i]) return nadds;
 
       size_t local_point_ids_array_size = 0;
       size_t num_local_point_ids;

src/grid_point_search.h

@@ -37,6 +37,7 @@ squareDistance(const double *a, const double *b) noexcept
 
 enum class PointSearchMethod
 {
+  undefined,
   full,
   nanoflann,
   kdtree,
@@ -51,7 +52,7 @@ struct GridPointSearch
   bool is_cyclic = false;
   bool is_reg2d = false;
   bool is_curve = false;
-  PointSearchMethod method = PointSearchMethod::full;
+  PointSearchMethod method = PointSearchMethod::nanoflann;
   size_t n = 0;
   size_t dims[2] = { 0 };
 
@@ -78,6 +79,7 @@ void gridSearchPointSmooth(GridPointSearch &gps, double plon, double plat, knnWe
 
 void gridPointSearchCreateReg2d(GridPointSearch &gps, bool xIsCyclic, size_t dims[2], const Varray<double> &lons, const Varray<double> &lats);
 void gridPointSearchCreate(GridPointSearch &gps, bool xIsCyclic, size_t dims[2], size_t n, const Varray<double> &lons, const Varray<double> &lats);
+void gridPointSearchCreate(GridPointSearch &gps, const Varray<double> &lons, const Varray<double> &lats, PointSearchMethod method = PointSearchMethod::spherepart);
 void gridPointSearchDelete(GridPointSearch &gps);
 size_t gridPointSearchNearest(GridPointSearch &gps, double lon, double lat, size_t *addr, double *dist);
 size_t gridPointSearchQnearest(GridPointSearch &gps, double lon, double lat, size_t nnn, size_t *adds, double *dist);

src/remap_point_search.cc

@@ -43,7 +43,7 @@ fill_src_add(bool is_cyclic, long nx, long ny, long ii, long jj, long k, size_t
   for (long j = j0; j <= jn; ++j)
     for (long i = i0; i <= in; ++i)
       {
-        long ix = i;
+        auto ix = i;
         if (is_cyclic && ix < 0) ix += nx;
         if (is_cyclic && ix >= nx) ix -= nx;
         if (ix >= 0 && ix < nx && j >= 0 && j < ny) psrc_add[num_add++] = j * nx + ix;

src/remaplib.cc

@@ -630,17 +630,17 @@ remapSearchInit(RemapMethod mapType, RemapSearch &search, RemapGrid &src_grid, R
   search.srcBins.nbins = remap_num_srch_bins;
   search.tgtBins.nbins = remap_num_srch_bins;
 
-  bool usePointsearch = mapType == RemapMethod::DISTWGT;
+  auto usePointsearch = (mapType == RemapMethod::DISTWGT);
   if (src_grid.type != RemapGridType::Reg2D && pointSearchMethod != PointSearchMethod::latbins)
     {
-      usePointsearch |= mapType == RemapMethod::BILINEAR;
-      usePointsearch |= mapType == RemapMethod::BICUBIC;
+      usePointsearch |= (mapType == RemapMethod::BILINEAR);
+      usePointsearch |= (mapType == RemapMethod::BICUBIC);
     }
 
-  bool useCellsearch = (mapType == RemapMethod::CONSERV)
+  auto useCellsearch = (mapType == RemapMethod::CONSERV)
                        && (cellSearchMethod == CellSearchMethod::spherepart || src_grid.type == RemapGridType::Reg2D);
 
-  double start = Options::cdoVerbose ? cdo_get_wtime() : 0;
+  auto start = Options::cdoVerbose ? cdo_get_wtime() : 0.0;
 
   if (usePointsearch)
     {