Remapstat: test implementation for unstructured target grids using yac_point_in_cell().

src/Remapstat.cc

@@ -33,6 +33,15 @@
 #include "verifygrid.h"
 
 //#define TESTIMPL 1
+//#define USE_YAC  1
+#ifdef TESTIMPL
+#ifdef USE_YAC
+extern "C"
+{
+#include "lib/yac/geometry.h"
+}
+#endif
+#endif
 
 constexpr double PI2 = M_PI * 2.0;
 
@@ -233,6 +242,66 @@ double calc_maxdist_rec2d(size_t i, size_t nlon, double plon, double plat, const
   return maxdist;
 }
 
+#ifdef USE_YAC
+static
+size_t find_points_yac(std::vector<char> &vmask, size_t cell_no, size_t ncorner, size_t nadds, std::vector<size_t> &adds, const Varray<double> &xvals1,
+                   const Varray<double> &yvals1, const Varray<double> &xbounds, const Varray<double> &ybounds)
+{
+#ifndef TESTIMPL
+  cdoAbort("Internal error: find_points() not implemented!");
+#endif
+
+  struct grid_cell cell;
+  cell.array_size = ncorner;
+  cell.num_corners = ncorner;
+  cell.edge_type = new enum yac_edge_type[ncorner];
+  cell.coordinates_x = new double[ncorner];
+  cell.coordinates_y = new double[ncorner];
+  cell.coordinates_xyz = new double[ncorner][3];
+
+  // enum yac_edge_type lonlat_circle_type[] = { LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE };
+
+  for (size_t k = 0; k < ncorner; ++k)
+    {
+      const auto lon = xbounds[cell_no * ncorner + k];
+      const auto lat = ybounds[cell_no * ncorner + k];
+      cell.edge_type[k] = GREAT_CIRCLE;
+      // cell.edge_type[k] = lonlat_circle_type[k + 1];
+      cell.coordinates_x[k] = lon;
+      cell.coordinates_y[k] = lat;
+      gcLLtoXYZ(lon, lat, cell.coordinates_xyz[k]);
+    }
+
+  double center_coordinates[3];
+  size_t nvalues = 0;
+  for (size_t k = 0; k < nadds; ++k)
+    {
+      const auto index1 = adds[k];
+      const auto lon = xvals1[index1];
+      const auto lat = yvals1[index1];
+
+      gcLLtoXYZ(lon, lat, center_coordinates);
+      
+      if (yac_point_in_cell(center_coordinates, cell))
+        {
+          if (vmask[index1] == 0)
+            {
+              adds[nvalues] = adds[k];
+              nvalues++;
+            }
+          if (vmask[index1] < 127) vmask[index1]++;
+        }
+    }
+
+  delete[] cell.edge_type;
+  delete[] cell.coordinates_x;
+  delete[] cell.coordinates_y;
+  delete[] cell.coordinates_xyz;
+
+  return nvalues;
+}
+#endif
+
 static
 size_t find_points(std::vector<char> &vmask, size_t cell_no, size_t ncorner, size_t nadds, std::vector<size_t> &adds, const Varray<double> &xvals1,
                    const Varray<double> &yvals1, const Varray<double> &xbounds, const Varray<double> &ybounds)
@@ -326,9 +395,12 @@ size_t find_points_rec2d(std::vector<char> &vmask, size_t i, size_t nlon2, size_
         if ((x >= xbounds2d[2*ix] && x < xbounds2d[2*ix + 1]) ||
             ((x-PI2) >= xbounds2d[2*ix] && (x-PI2) < xbounds2d[2*ix + 1]))
           {
+            if (vmask[index1] == 0)
+              {
+                adds[nvalues] = index1;
+                nvalues++;
+              }
             if (vmask[index1] < 127) vmask[index1]++;
-            adds[nvalues] = index1;
-            nvalues++;
           }
     }
 
@@ -426,7 +498,11 @@ Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
 
       const auto nvalues = grid2_is_reg2d
         ? find_points_rec2d(vmask, i, nlon2, nadds, adds, xvals1, yvals1, xbounds2d, ybounds2d)
+#ifdef USE_YAC
+        : find_points_yac(vmask, i, nv, nadds, adds, xvals1, yvals1, xbounds, ybounds);
+#else
         : find_points(vmask, i, nv, nadds, adds, xvals1, yvals1, xbounds, ybounds);
+#endif
 
       if (nvalues)
         {