Added function gridSearchPointHealpix().

src/remap_point_search.cc

@@ -43,8 +43,41 @@ fill_src_add(bool isCyclic, long nx, long ny, long ii, long jj, long k, size_t *
 }
 
 static void
-store_distance(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights, size_t nx, size_t num_add,
-               size_t *psrc_add)
+store_distance_healpix(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights, size_t num_add, size_t *psrc_add)
+{
+  double xyz[3], query_pt[3];
+  gcLLtoXYZ(plon, plat, query_pt);
+  auto sqrSearchRadius = cdo::sqr(gps.searchRadius);
+
+  for (size_t na = 0; na < num_add; ++na)
+    {
+      auto nadd = psrc_add[na];
+      /*
+      auto iy = nadd / nx;
+      auto ix = nadd - iy * nx;
+
+      xyz[0] = coslat[iy] * coslon[ix];
+      xyz[1] = coslat[iy] * sinlon[ix];
+      xyz[2] = sinlat[iy];
+      */
+      xyz[0] = query_pt[0];
+      xyz[1] = query_pt[1];
+      xyz[2] = query_pt[2];
+      // Find distance to this point
+      double sqrDist = (float) squareDistance(query_pt, xyz);
+      if (sqrDist <= sqrSearchRadius)
+        {
+          // Store the address and distance if this is one of the smallest so far
+          knnWeights.storeDistance(nadd, std::sqrt(sqrDist));
+        }
+    }
+
+  knnWeights.checkDistance();
+}
+
+static void
+store_distance_reg2d(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights, size_t nx, size_t num_add,
+                     size_t *psrc_add)
 {
   const auto &coslon = gps.coslon;
   const auto &sinlon = gps.sinlon;
@@ -53,19 +86,19 @@ store_distance(GridPointSearch &gps, double plon, double plat, knnWeightsType &k
 
   double xyz[3], query_pt[3];
   gcLLtoXYZ(plon, plat, query_pt);
-  const auto sqrSearchRadius = cdo::sqr(gps.searchRadius);
+  auto sqrSearchRadius = cdo::sqr(gps.searchRadius);
 
   for (size_t na = 0; na < num_add; ++na)
     {
-      const auto nadd = psrc_add[na];
-      const auto iy = nadd / nx;
-      const auto ix = nadd - iy * nx;
+      auto nadd = psrc_add[na];
+      auto iy = nadd / nx;
+      auto ix = nadd - iy * nx;
 
       xyz[0] = coslat[iy] * coslon[ix];
       xyz[1] = coslat[iy] * sinlon[ix];
       xyz[2] = sinlat[iy];
       // Find distance to this point
-      const double sqrDist = (float) squareDistance(query_pt, xyz);
+      double sqrDist = (float) squareDistance(query_pt, xyz);
       if (sqrDist <= sqrSearchRadius)
         {
           // Store the address and distance if this is one of the smallest so far
@@ -76,6 +109,32 @@ store_distance(GridPointSearch &gps, double plon, double plat, knnWeightsType &k
   knnWeights.checkDistance();
 }
 
+static void
+gridSearchPointHealpix(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights)
+{
+  /*
+    Input variables:
+
+      plat : latitude  of the search point
+      plon : longitude of the search point
+
+    Output variables:
+
+      knnWeights.m_addr[numNeighbors] :  address of each of the closest points
+      knnWeights.m_dist[numNeighbors] : distance to each of the closest points
+  */
+  auto numNeighbors = knnWeights.maxNeighbors();
+
+  // Initialize distance and address arrays
+  knnWeights.initAddr();
+  knnWeights.initDist();
+
+  auto index = hp_lonlat_to_healpixl(gps.order, gps.nside, plon, plat);
+
+  size_t src_add = index;
+  store_distance_healpix(gps, plon, plat, knnWeights, 1, &src_add);
+}
+
 // This routine finds the closest numNeighbor points to a search point and computes a distance to each of the neighbors
 static void
 gridSearchPointReg2d(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights)
@@ -91,7 +150,7 @@ gridSearchPointReg2d(GridPointSearch &gps, double plon, double plat, knnWeightsT
       knnWeights.m_addr[numNeighbors] :  address of each of the closest points
       knnWeights.m_dist[numNeighbors] : distance to each of the closest points
   */
-  const auto numNeighbors = knnWeights.maxNeighbors();
+  auto numNeighbors = knnWeights.maxNeighbors();
   auto &nbr_add = knnWeights.m_addr;
   auto &nbr_dist = knnWeights.m_dist;
 
@@ -130,9 +189,9 @@ gridSearchPointReg2d(GridPointSearch &gps, double plon, double plat, knnWeightsT
           psrc_add = src_add_tmp.data();
         }
 
-      const auto num_add = fill_src_add(gps.isCyclic, nx, ny, ii, jj, k, psrc_add);
+      auto num_add = fill_src_add(gps.isCyclic, nx, ny, ii, jj, k, psrc_add);
 
-      store_distance(gps, plon, plat, knnWeights, nx, num_add, psrc_add);
+      store_distance_reg2d(gps, plon, plat, knnWeights, nx, num_add, psrc_add);
     }
   else if (gps.extrapolate)
     {
@@ -218,7 +277,7 @@ grid_search_point_smooth(GridPointSearch &gps, double plon, double plat, knnWeig
       knnWeights.m_dist[numNeighbors] : distance to each of the closest points
   */
   auto numNeighbors = knnWeights.maxNeighbors();
-  const auto checkDistance = (numNeighbors <= 32);
+  auto checkDistance = (numNeighbors <= 32);
 
   // check some more points if distance is the same use the smaller index (nadd)
   auto ndist = checkDistance ? ((numNeighbors > 8) ? numNeighbors + 8 : numNeighbors * 2) : numNeighbors;
@@ -259,7 +318,9 @@ grid_search_point_smooth(GridPointSearch &gps, double plon, double plat, knnWeig
 void
 remap_search_points(RemapSearch &rsearch, const LonLatPoint &llp, knnWeightsType &knnWeights)
 {
-  if (rsearch.srcGrid->type == RemapGridType::Reg2D)
+  if (rsearch.srcGrid->type == RemapGridType::HealPix)
+    gridSearchPointHealpix(rsearch.gps, llp.lon, llp.lat, knnWeights);
+  else if (rsearch.srcGrid->type == RemapGridType::Reg2D)
     gridSearchPointReg2d(rsearch.gps, llp.lon, llp.lat, knnWeights);
   else
     grid_search_point(rsearch.gps, llp.lon, llp.lat, knnWeights);
@@ -290,8 +351,8 @@ gridSearchSquareCurv2d(GridPointSearch &gps, RemapGrid *rgrid, size_t (&src_add)
   size_t nadds = grid_point_search_nearest(gps, plon, plat, &addr, &dist);
   if (nadds > 0)
     {
-      const auto nx = rgrid->dims[0];
-      const auto ny = rgrid->dims[1];
+      auto nx = rgrid->dims[0];
+      auto ny = rgrid->dims[1];
 
       for (int k = 0; k < 4; ++k)
         {