Changed interface of gridsearch_create_nn().

src/Fillmiss.cc

@@ -346,10 +346,12 @@ void setmisstodis(field_type *field1, field_type *field2, int num_neighbors)
 
   if ( nmiss )
     {
+      bool xIsCyclic = false;
+      size_t dims[2] = {nvals, 0};
       if ( num_neighbors == 1 )
-        gs = gridsearch_create_nn(nvals, lons, lats);
+        gs = gridsearch_create_nn(xIsCyclic, dims, nvals, lons, lats);
       else
-        gs = gridsearch_create(nvals, lons, lats);
+        gs = gridsearch_create(xIsCyclic, dims, nvals, lons, lats);
     }
   
   finish = clock();

src/Samplegridicon.cc

@@ -299,7 +299,9 @@ static
 void compute_child_from_bounds(cellindex_type *cellindex2, long ncells2, double *grid_center_lon2, double *grid_center_lat2, double *grid_corner_lon2,
                                double *grid_corner_lat2, long ncells1, double *grid_center_lon1, double *grid_center_lat1)
 {
-  struct gridsearch *gs = gridsearch_create(ncells1, grid_center_lon1, grid_center_lat1);
+  bool xIsCyclic = false;
+  size_t dims[2] = {ncells1, 0};
+  struct gridsearch *gs = gridsearch_create(xIsCyclic, dims, ncells1, grid_center_lon1, grid_center_lat1);
   size_t nbr_add[MAX_SEARCH];  // source address at nearest neighbors
   double nbr_dist[MAX_SEARCH]; // angular distance four nearest neighbors
 

src/Smooth.cc

@@ -133,10 +133,12 @@ void smooth(int gridID, double missval, const double *restrict array1, double *r
 
   struct gridsearch *gs = NULL;
 
+  bool xIsCyclic = false;
+  size_t dims[2] = {gridsize, 0};
   if ( num_neighbors == 1 )
-    gs = gridsearch_create_nn(gridsize, xvals, yvals);
+    gs = gridsearch_create_nn(xIsCyclic, dims, gridsize, xvals, yvals);
   else
-    gs = gridsearch_create(gridsize, xvals, yvals);
+    gs = gridsearch_create(xIsCyclic, dims, gridsize, xvals, yvals);
 
   gs->search_radius = spoint.radius;
 

src/grid_search.cc

@@ -134,11 +134,11 @@ void gridsearch_extrapolate(struct gridsearch *gs)
 }
 
 
-struct gridsearch *gridsearch_create_reg2d(bool is_cyclic, size_t dims[2], const double *restrict lons, const double *restrict lats)
+struct gridsearch *gridsearch_create_reg2d(bool xIsCyclic, size_t dims[2], const double *restrict lons, const double *restrict lats)
 {
   struct gridsearch *gs = (struct gridsearch *) Calloc(1, sizeof(struct gridsearch));
 
-  gs->is_cyclic = is_cyclic;
+  gs->is_cyclic = xIsCyclic;
   gs->is_reg2d = true;
   gs->dims[0] = dims[0];
   gs->dims[1] = dims[1];
@@ -146,7 +146,7 @@ struct gridsearch *gridsearch_create_reg2d(bool is_cyclic, size_t dims[2], const
   size_t ny = dims[0];
 
   size_t nxm = nx;
-  if ( is_cyclic ) nxm++;
+  if ( xIsCyclic ) nxm++;
 
   double *reg2d_center_lon = (double *) Malloc(nxm*sizeof(double));
   double *reg2d_center_lat = (double *) Malloc(ny*sizeof(double));
@@ -490,10 +490,15 @@ void gridsearch_bound_poly(struct gridsearch *gs, size_t dims[2], size_t n, cons
 }
 
 
-struct gridsearch *gridsearch_create(size_t n, const double *restrict lons, const double *restrict lats)
+struct gridsearch *gridsearch_create(bool xIsCyclic, size_t dims[2], size_t n, const double *restrict lons, const double *restrict lats)
 {
   struct gridsearch *gs = (struct gridsearch *) Calloc(1, sizeof(struct gridsearch));
 
+  gs->is_cyclic = xIsCyclic;
+  gs->is_curve = n!=1 && n==dims[0]*dims[1];
+  gs->dims[0] = dims[0];
+  gs->dims[1] = dims[1];
+
   gs->method_nn = gridsearch_method_nn;
   gs->n = n;
   if ( n == 0 ) return gs;
@@ -507,10 +512,15 @@ struct gridsearch *gridsearch_create(size_t n, const double *restrict lons, cons
 }
 
 
-struct gridsearch *gridsearch_create_nn(size_t n, const double *restrict lons, const double *restrict lats)
+struct gridsearch *gridsearch_create_nn(bool xIsCyclic, size_t dims[2], size_t n, const double *restrict lons, const double *restrict lats)
 {
   struct gridsearch *gs = (struct gridsearch *) Calloc(1, sizeof(struct gridsearch));
 
+  gs->is_cyclic = xIsCyclic;
+  gs->is_curve = n!=1 && n==dims[0]*dims[1];
+  gs->dims[0] = dims[0];
+  gs->dims[1] = dims[1];
+
   gs->method_nn = gridsearch_method_nn;
   gs->n = n;
   if ( n == 0 ) return gs;
@@ -629,33 +639,39 @@ size_t gs_nearest_nanoflann(void *search_container, double lon, double lat, doub
   //printf("%zu %g\n", ret_index, out_dist_sqr);
 
 #ifdef  TEST_QUAD
-  if ( ret_index != GS_NOT_FOUND )
+  if ( !gs->extrapolate && gs->is_curve )
     {
-      size_t nx = 450;
-      size_t ny = 250;
-      size_t adds[4];
-      double lons[4];
-      double lats[4];
-      bool is_cyclic = false;
-      for ( unsigned k = 0; k < 4; ++k )
+      if ( ret_index != GS_NOT_FOUND )
         {
-          /* Determine neighbor addresses */
-          size_t j = ret_index/nx;
-          size_t i = ret_index - j*nx;
-          if ( k == 1 || k == 3 )  i = (i > 0) ? i - 1 : (is_cyclic) ? nx-1 : 0;
-          if ( k == 2 || k == 3 )  j = (j > 0) ? j - 1 : 0;
-
-          if ( point_in_quad(is_cyclic, nx, ny, i, j, adds, lons, lats,
-                             lon, lat, gs->plons, gs->plats) )
+          size_t nx = gs->dims[0];
+          size_t ny = gs->dims[1];
+          size_t adds[4];
+          double lons[4];
+          double lats[4];
+          bool is_cyclic = gs->is_cyclic;
+          for ( unsigned k = 0; k < 4; ++k )
             {
-              index = ret_index;
-              break;
+              /* Determine neighbor addresses */
+              size_t j = ret_index/nx;
+              size_t i = ret_index - j*nx;
+              if ( k == 1 || k == 3 )  i = (i > 0) ? i - 1 : (is_cyclic) ? nx-1 : 0;
+              if ( k == 2 || k == 3 )  j = (j > 0) ? j - 1 : 0;
+
+              if ( point_in_quad(is_cyclic, nx, ny, i, j, adds, lons, lats,
+                                 lon, lat, gs->plons, gs->plats) )
+                {
+                  index = ret_index;
+                  break;
+                }
             }
         }
     }
+  else
 #else
-  index = ret_index;
-  *prange = out_dist_sqr; 
+    {
+      index = ret_index;
+      *prange = out_dist_sqr;
+    }
 #endif
 
   //float frange = range;

src/grid_search.h

@@ -12,6 +12,7 @@ struct gridsearch {
   bool extrapolate;
   bool is_cyclic;
   bool is_reg2d;
+  bool is_curve;
   GridsearchMethod method_nn;
   size_t n;
   size_t dims[2];
@@ -44,9 +45,9 @@ struct gsknn *gridsearch_knn_new(size_t size);
 void gridsearch_knn_delete(struct gsknn *knn);
 size_t gridsearch_knn(struct gridsearch *gs, struct gsknn *knn, double plon, double plat);
 
-struct gridsearch *gridsearch_create_reg2d(bool is_cyclic, size_t dims[2], const double *restrict lons, const double *restrict lats);
-struct gridsearch *gridsearch_create(size_t n, const double *restrict lons, const double *restrict lats);
-struct gridsearch *gridsearch_create_nn(size_t n, const double *restrict lons, const double *restrict lats);
+struct gridsearch *gridsearch_create_reg2d(bool xIsCyclic, size_t dims[2], const double *restrict lons, const double *restrict lats);
+struct gridsearch *gridsearch_create(bool xIsCyclic, size_t dims[2], size_t n, const double *restrict lons, const double *restrict lats);
+struct gridsearch *gridsearch_create_nn(bool xIsCyclic, size_t dims[2], size_t n, const double *restrict lons, const double *restrict lats);
 void gridsearch_delete(struct gridsearch *gs);
 size_t gridsearch_nearest(struct gridsearch *gs, double lon, double lat, double *range);
 size_t gridsearch_qnearest(struct gridsearch *gs, double lon, double lat, double *prange, size_t nnn, size_t *adds, double *dist);

src/remap_bilinear_scrip.cc

@@ -371,9 +371,11 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
   progressInit();
 
 #ifdef  TEST_KDTREE
+  bool xIsCyclic = false;
+  size_t dims[2] = {src_grid->size, 0};
   struct gridsearch *gs = NULL;
   if (remap_grid_type != REMAP_GRID_TYPE_REG2D )
-    gs = gridsearch_create_nn(src_grid->size, src_grid->cell_center_lon, src_grid->cell_center_lat);
+    gs = gridsearch_create_nn(xIsCyclic, dims, src_grid->size, src_grid->cell_center_lon, src_grid->cell_center_lat);
 #endif
 
   size_t tgt_grid_size = tgt_grid->size;

src/remap_distwgt.cc

@@ -366,13 +366,15 @@ void remap_distwgt_weights(size_t num_neighbors, remapgrid_t *src_grid, remapgri
   double start = cdoVerbose ? omp_get_wtime() : 0;
 #endif
 
+  bool xIsCyclic = src_grid->is_cyclic;
+  size_t *dims = src_grid->dims;
   struct gridsearch *gs = NULL;
   if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
-    gs = gridsearch_create_reg2d(src_grid->is_cyclic, src_grid->dims, src_grid->reg2d_center_lon, src_grid->reg2d_center_lat);
+    gs = gridsearch_create_reg2d(xIsCyclic, dims, src_grid->reg2d_center_lon, src_grid->reg2d_center_lat);
   else if ( num_neighbors == 1 )
-    gs = gridsearch_create_nn(src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
+    gs = gridsearch_create_nn(xIsCyclic, dims, src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
   else
-    gs = gridsearch_create(src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
+    gs = gridsearch_create(xIsCyclic, dims, src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
 
   if ( src_grid->lextrapolate ) gridsearch_extrapolate(gs);
   //else
@@ -475,13 +477,15 @@ void remap_distwgt(size_t num_neighbors, remapgrid_t *src_grid, remapgrid_t *tgt
   double start = cdoVerbose ? omp_get_wtime() : 0;
 #endif
 
+  bool xIsCyclic = src_grid->is_cyclic;
+  size_t *dims = src_grid->dims;
   struct gridsearch *gs = NULL;
   if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
-    gs = gridsearch_create_reg2d(src_grid->is_cyclic, src_grid->dims, src_grid->reg2d_center_lon, src_grid->reg2d_center_lat);
+    gs = gridsearch_create_reg2d(xIsCyclic, dims, src_grid->reg2d_center_lon, src_grid->reg2d_center_lat);
   else if ( num_neighbors == 1 )
-    gs = gridsearch_create_nn(src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
+    gs = gridsearch_create_nn(xIsCyclic, dims, src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
   else
-    gs = gridsearch_create(src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
+    gs = gridsearch_create(xIsCyclic, dims, src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
 
   if ( src_grid->lextrapolate ) gridsearch_extrapolate(gs);
 
@@ -611,13 +615,15 @@ void intgriddis(field_type *field1, field_type *field2, size_t num_neighbors)
   double start = cdoVerbose ? omp_get_wtime() : 0;
 #endif
 
+  bool xIsCyclic = gridIsCircular(gridID1);
+  size_t dims[2] = {src_grid_size, 0};
   struct gridsearch *gs = NULL;
   // if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
-  //  gs = gridsearch_create_reg2d(src_grid->is_cyclic, src_grid->dims, src_grid->reg2d_center_lon, src_grid->reg2d_center_lat);
+  //  gs = gridsearch_create_reg2d(xIsCyclic, dims, src_grid->reg2d_center_lon, src_grid->reg2d_center_lat);
   if ( num_neighbors == 1 )
-    gs = gridsearch_create_nn(src_grid_size, src_cell_center_lon, src_cell_center_lat);
+    gs = gridsearch_create_nn(xIsCyclic, dims, src_grid_size, src_cell_center_lon, src_cell_center_lat);
   else
-    gs = gridsearch_create(src_grid_size, src_cell_center_lon, src_cell_center_lat);
+    gs = gridsearch_create(xIsCyclic, dims, src_grid_size, src_cell_center_lon, src_cell_center_lat);
 
   // if ( src_grid->lextrapolate ) gridsearch_extrapolate(gs);
   // gridsearch_extrapolate(gs);