Added gridsearch_extrapolate().

src/Remap.cc

@@ -117,7 +117,7 @@ int maptype2operfunc(int map_type, int submap_type, int num_neighbors, int remap
 } 
 
 static
-void print_remap_info(int operfunc, int remap_genweights, remapgrid_t *src_grid, remapgrid_t *tgt_grid, size_t nmiss)
+void print_remap_info(int operfunc, bool remap_genweights, remapgrid_t *src_grid, remapgrid_t *tgt_grid, size_t nmiss)
 {
   char line[256], tmpstr[256];
 
@@ -218,7 +218,7 @@ static bool lextrapolate = false;
 int max_remaps = -1;
 int sort_mode = HEAP_SORT;
 double remap_frac_min = 0;
-int remap_genweights = TRUE;
+bool REMAP_genweights = true;
 
 static
 void get_remap_env(void)
@@ -417,19 +417,15 @@ void get_remap_env(void)
       if ( *envstr )
 	{
 	  if ( memcmp(envstr, "ON", 2) == 0 || memcmp(envstr, "on", 2) == 0 )
-	    {
-	      remap_genweights = TRUE;
-	    }
+            REMAP_genweights = true;
 	  else if ( memcmp(envstr, "OFF", 3) == 0 || memcmp(envstr, "off", 3) == 0 )
-	    {
-	      remap_genweights = FALSE;
-	    }
+            REMAP_genweights = false;
 	  else
 	    cdoWarning("Environment variable CDO_REMAP_GENWEIGHTS has wrong value!");
 
 	  if ( cdoVerbose )
 	    {
-	      if ( remap_genweights == TRUE )
+	      if ( REMAP_genweights )
 		cdoPrint("Generation of weights enabled!");
 	      else
 		cdoPrint("Generation of weights disabled!");
@@ -775,6 +771,8 @@ void sort_remap_add(remapvars_t *remapvars)
 
 void *Remap(void *argument)
 {
+  bool remap_genweights = REMAP_genweights;
+  bool need_gradiants = false;
   int streamID2 = -1;
   int nrecs;
   int varID, levelID;
@@ -788,7 +786,6 @@ void *Remap(void *argument)
   int map_type = -1;
   int submap_type = SUBMAP_TYPE_NONE;
   int num_neighbors = 0;
-  int need_gradiants = FALSE;
   char varname[CDI_MAX_NAME];
   double missval;
   remap_t *remaps = NULL;
@@ -901,8 +898,7 @@ void *Remap(void *argument)
 	}
     }
 
-  if ( lwrite_remap || lremapxxx )
-    remap_genweights = TRUE;
+  if ( lwrite_remap || lremapxxx ) remap_genweights = true;
 
   if ( lremapxxx )
     {
@@ -963,22 +959,22 @@ void *Remap(void *argument)
       get_map_type(operfunc, &map_type, &submap_type, &num_neighbors, &remap_order);
     }
 
-  if ( remap_genweights == FALSE &&
+  if ( !remap_genweights &&
        map_type != MAP_TYPE_BILINEAR && map_type != MAP_TYPE_BICUBIC &&
        map_type != MAP_TYPE_DISTWGT  && map_type != MAP_TYPE_CONSERV_YAC )
-    remap_genweights = TRUE;
+    remap_genweights = true;
 
-  remap_set_int(REMAP_GENWEIGHTS, remap_genweights);
+  remap_set_int(REMAP_GENWEIGHTS, (int)remap_genweights);
 
   if ( map_type == MAP_TYPE_CONSERV || map_type == MAP_TYPE_CONSERV_YAC ) norm_opt = get_norm_opt();
 
   size_t grid1sizemax = vlistGridsizeMax(vlistID1);
 
-  if ( map_type == MAP_TYPE_BICUBIC ) need_gradiants = TRUE;
+  if ( map_type == MAP_TYPE_BICUBIC ) need_gradiants = true;
   if ( map_type == MAP_TYPE_CONSERV && remap_order == 2 )
     {
       if ( cdoVerbose ) cdoPrint("Second order remapping");
-      need_gradiants = TRUE;
+      need_gradiants = true;
     }
   else
     remap_order = 1;
@@ -1043,13 +1039,6 @@ void *Remap(void *argument)
 		  grid1sizemax = gridsize_new;
 		  array1 = (double*) Realloc(array1, grid1sizemax*sizeof(double));
 		  imask  = (int*) Realloc(imask, grid1sizemax*sizeof(int));
-
-		  if ( need_gradiants )
-		    {
-		      grad1_lat    = (double*) Realloc(grad1_lat, grid1sizemax*sizeof(double));
-		      grad1_lon    = (double*) Realloc(grad1_lon, grid1sizemax*sizeof(double));
-		      grad1_latlon = (double*) Realloc(grad1_latlon, grid1sizemax*sizeof(double));
-		    }
 		}
 	      
 	      for ( long j = ny-1; j >= 0; j-- )
@@ -1106,9 +1095,7 @@ void *Remap(void *argument)
 		  if ( gridIsCircular(gridID1) && !lextrapolate ) remap_extrapolate = true;
 		  remaps[r].src_grid.non_global = false;
 		  if ( map_type == MAP_TYPE_DISTWGT && !remap_extrapolate && gridInqSize(gridID1) > 1 && !is_global_grid(gridID1) )
-		    {
-		      remaps[r].src_grid.non_global = true;
-		    }
+                    remaps[r].src_grid.non_global = true;
 		  /*
 		    remaps[r].src_grid.luse_cell_area = FALSE;
 		    remaps[r].tgt_grid.luse_cell_area = FALSE;
@@ -1205,7 +1192,7 @@ void *Remap(void *argument)
 	      if ( need_gradiants )
 		{
 		  if ( remaps[r].src_grid.rank != 2 && remap_order == 2 )
-		    cdoAbort("Second order remapping is not only available for unstructured grids!");
+		    cdoAbort("Second order remapping is not available for unstructured grids!");
 
 		  remap_gradients(remaps[r].src_grid, array1, grad1_lat, grad1_lon, grad1_latlon);
 		}

src/cdo.cc

@@ -1258,7 +1258,7 @@ int parse_options_long(int argc, char *argv[])
               int intarg = parameter2int(CDO_optarg);
               if ( intarg != 0 && intarg != 1 )
                 cdoAbort("Unsupported value for option --remap_genweights %d [0/1]", intarg);
-              remap_genweights = intarg;
+              REMAP_genweights = intarg;
             }
           else if ( lsortname )
             {

src/grid_search.cc

@@ -77,10 +77,17 @@ void gridsearch_set_method(const char *methodstr)
 }
 
 
+void gridsearch_extrapolate(struct gridsearch *gs)
+{
+  gs->extrapolate = true;
+}
+
+
 struct gridsearch *gridsearch_create_reg2d(bool lcyclic, size_t nx, size_t ny, const double *restrict lons, const double *restrict lats)
 {
   struct gridsearch *gs = (struct gridsearch *) Calloc(1, sizeof(struct gridsearch));
 
+  gs->reg2d = true;
   gs->nx = nx;
   gs->ny = ny;
 

src/grid_search.h

@@ -26,6 +26,8 @@ struct gsNear {
 };
 
 struct gridsearch {
+  bool extrapolate;
+  bool reg2d;
   int method_nn;
   size_t n;
   size_t nx, ny;
@@ -62,5 +64,6 @@ struct gridsearch *gridsearch_create_nn(size_t n, const double *restrict lons, c
 void gridsearch_delete(struct gridsearch *gs);
 size_t gridsearch_nearest(struct gridsearch *gs, double lon, double lat, double *range);
 struct pqueue *gridsearch_qnearest(struct gridsearch *gs, double lon, double lat, double *prange, size_t nnn);
+void gridsearch_extrapolate(struct gridsearch *gs);
 
 #endif

src/remap_distwgt.cc

@@ -127,14 +127,13 @@ void grid_search_nbr_reg2d(struct gridsearch *gs, size_t num_neighbors, remapgri
     double plat,         ! latitude  of the search point
     double plon,         ! longitude of the search point
   */
-  size_t n, nadd;
+  size_t n;
   size_t ii, jj;
   long i, j, ix;
   size_t src_add[MAX_SEARCH_CELLS];
   size_t *src_add_tmp = NULL;
   size_t *psrc_add = src_add;
   size_t num_add = 0;
-  double distance;   //  Angular distance
   double cos_search_radius = cos(gs->search_radius);
   double coslat_dst = cos(plat);  // cos(lat)  of the search point
   double coslon_dst = cos(plon);  // cos(lon)  of the search point
@@ -204,11 +203,12 @@ void grid_search_nbr_reg2d(struct gridsearch *gs, size_t num_neighbors, remapgri
   if ( lfound )
     {
       size_t ix, iy;
+      size_t nadd;
+      double distance;   //  Angular distance
 
       for ( size_t na = 0; na < num_add; ++na )
 	{
 	  nadd = psrc_add[na];
-
 	  iy = nadd/nx;
 	  ix = nadd - iy*nx;
 
@@ -232,7 +232,7 @@ void grid_search_nbr_reg2d(struct gridsearch *gs, size_t num_neighbors, remapgri
 
       if ( src_add_tmp ) Free(src_add_tmp);
     }
-  else if ( src_grid->lextrapolate )
+  else if ( gs->extrapolate )
     {
       int search_result;
 
@@ -397,6 +397,8 @@ void remap_distwgt_weights(size_t num_neighbors, remapgrid_t *src_grid, remapgri
   else
     gs = gridsearch_create(src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
 
+  if ( src_grid->lextrapolate ) gridsearch_extrapolate(gs);
+
 #if defined(_OPENMP)
   if ( cdoVerbose ) printf("gridsearch created: %.2f seconds\n", omp_get_wtime()-start);
   if ( cdoVerbose ) start = omp_get_wtime();
@@ -506,6 +508,8 @@ void remap_distwgt(size_t num_neighbors, remapgrid_t *src_grid, remapgrid_t *tgt
   else
     gs = gridsearch_create(src_grid_size, src_grid->cell_center_lon, src_grid->cell_center_lat);
 
+  if ( src_grid->lextrapolate ) gridsearch_extrapolate(gs);
+
 #if defined(_OPENMP)
   if ( cdoVerbose ) printf("gridsearch created: %.2f seconds\n", omp_get_wtime()-start);
   if ( cdoVerbose ) start = omp_get_wtime();

src/remap_search_reg2d.cc

@@ -6,20 +6,13 @@ int grid_search_reg2d_nn(size_t nx, size_t ny, size_t *restrict nbr_add, double
 			 const double *restrict src_center_lat, const double *restrict src_center_lon)
 {
   int search_result = 0;
-  size_t n, srch_add;
-  size_t i;
-  size_t ii, jj;
-  size_t jjskip;
-  double coslat, sinlat;
-  double dist_min, distance; /* For computing dist-weighted avg */
-  double *sincoslon;
   double coslat_dst = cos(plat);
   double sinlat_dst = sin(plat);
   double coslon_dst = cos(plon);
   double sinlon_dst = sin(plon);
 
-  dist_min = BIGNUM;
-  for ( n = 0; n < 4; ++n ) nbr_dist[n] = BIGNUM;  
+  double dist_min = BIGNUM;
+  for ( unsigned n = 0; n < 4; ++n ) nbr_dist[n] = BIGNUM;  
 
   size_t jjf = 0, jjl = ny-1;
   if ( plon >= src_center_lon[0] && plon <= src_center_lon[nx-1] )
@@ -40,33 +33,31 @@ int grid_search_reg2d_nn(size_t nx, size_t ny, size_t *restrict nbr_add, double
 	}
     }
 
-  sincoslon = (double*) Malloc(nx*sizeof(double));
+  double *sincoslon = (double*) Malloc(nx*sizeof(double));
 
-  for ( ii = 0; ii < nx; ++ii )
+  for ( size_t ii = 0; ii < nx; ++ii )
     sincoslon[ii] = coslon_dst*cos(src_center_lon[ii]) + sinlon_dst*sin(src_center_lon[ii]);
 
-  for ( jj = jjf; jj <= jjl; ++jj )
+  for ( size_t jj = jjf; jj <= jjl; ++jj )
     {
-      coslat = coslat_dst*cos(src_center_lat[jj]);
-      sinlat = sinlat_dst*sin(src_center_lat[jj]);
+      double coslat = coslat_dst*cos(src_center_lat[jj]);
+      double sinlat = sinlat_dst*sin(src_center_lat[jj]);
 
-      jjskip = jj > 1 && jj < (ny-2);
+      size_t jjskip = jj > 1 && jj < (ny-2);
 
-      for ( ii = 0; ii < nx; ++ii )
+      for ( size_t ii = 0; ii < nx; ++ii )
 	{
 	  if ( jjskip && ii > 1 && ii < (nx-2) ) continue;
 
-	  srch_add = jj*nx + ii;
-
-	  distance = acos(coslat*sincoslon[ii] + sinlat);
-
+	  double distance = acos(coslat*sincoslon[ii] + sinlat);
 	  if ( distance < dist_min )
 	    {
-	      for ( n = 0; n < 4; ++n )
+              size_t srch_add = jj*nx + ii;
+	      for ( unsigned n = 0; n < 4; ++n )
 		{
 		  if ( distance < nbr_dist[n] )
 		    {
-		      for ( i = 3; i > n; --i )
+		      for ( unsigned i = 3; i > n; --i )
 			{
 			  nbr_add [i] = nbr_add [i-1];
 			  nbr_dist[i] = nbr_dist[i-1];
@@ -84,12 +75,12 @@ int grid_search_reg2d_nn(size_t nx, size_t ny, size_t *restrict nbr_add, double
 
   Free(sincoslon);
 
-  for ( n = 0; n < 4; ++n ) nbr_dist[n] = ONE/(nbr_dist[n] + TINY);
-  distance = 0.0;
-  for ( n = 0; n < 4; ++n ) distance += nbr_dist[n];
-  for ( n = 0; n < 4; ++n ) nbr_dist[n] /= distance;
+  for ( unsigned n = 0; n < 4; ++n ) nbr_dist[n] = ONE/(nbr_dist[n] + TINY);
+  double distance = 0.0;
+  for ( unsigned n = 0; n < 4; ++n ) distance += nbr_dist[n];
+  for ( unsigned n = 0; n < 4; ++n ) nbr_dist[n] /= distance;
 
-  return (search_result);
+  return search_result;
 }
 
 
@@ -114,12 +105,9 @@ int grid_search_reg2d(remapgrid_t *src_grid, size_t *restrict src_add, double *r
     double src_center_lat[]        ! latitude  of each src grid center 
     double src_center_lon[]        ! longitude of each src grid center
   */
-  /*  Local variables */
   int search_result = 0;
-  long n;
-  size_t ii, iix, jj;
 
-  for ( n = 0; n < 4; ++n ) src_add[n] = 0;
+  for ( unsigned n = 0; n < 4; ++n ) src_add[n] = 0;
 
   size_t nx = src_grid_dims[0];
   size_t ny = src_grid_dims[1];
@@ -130,11 +118,12 @@ int grid_search_reg2d(remapgrid_t *src_grid, size_t *restrict src_add, double *r
   if ( /*plon < 0   &&*/ plon < src_center_lon[0]     ) plon += PI2;
   if ( /*plon > PI2 &&*/ plon > src_center_lon[nxm-1] ) plon -= PI2;
 
+  size_t ii, jj;
   int lfound = rect_grid_search(&ii, &jj, plon, plat, nxm, ny, src_center_lon, src_center_lat);
 
   if ( lfound )
     {
-      iix = ii;
+      size_t iix = ii;
       if ( src_grid->is_cyclic && iix == (nxm-1) ) iix = 0;
       src_add[0] = (jj-1)*nx+(ii-1);
       src_add[1] = (jj-1)*nx+(iix);
@@ -175,5 +164,5 @@ int grid_search_reg2d(remapgrid_t *src_grid, size_t *restrict src_add, double *r
   */
   search_result = grid_search_reg2d_nn(nx, ny, src_add, src_lats, plat, plon, src_center_lat, src_center_lon);
 
-  return (search_result);
+  return search_result;
 }  /* grid_search_reg2d */

src/remaplib.cc

@@ -726,14 +726,11 @@ void cell_bounding_boxes(remapgrid_t *grid, int remap_grid_basis)
 	}
       else /* full grid search */
 	{
-	  size_t gridsize;
-	  size_t i, i4;
-	  
-	  gridsize = grid->size;
-  
 	  if ( cdoVerbose ) cdoPrint("Grid: bounds missing -> full grid search!");
 
-	  for ( i = 0; i < gridsize; ++i )
+	  size_t gridsize = grid->size;
+	  size_t i4;
+	  for ( size_t i = 0; i < gridsize; ++i )
 	    {
 	      i4 = i<<2;
 	      grid->cell_bound_box[i4  ] = RESTR_SCALE(-PIH);
@@ -794,10 +791,7 @@ void remap_grids_init(int map_type, bool lextrapolate, int gridID1, remapgrid_t
       if ( map_type == MAP_TYPE_BILINEAR && src_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
     }
 
-  if ( lextrapolate )
-    src_grid->lextrapolate = true;
-  else
-    src_grid->lextrapolate = false;
+  src_grid->lextrapolate = lextrapolate;
 
   if ( map_type == MAP_TYPE_CONSERV || map_type == MAP_TYPE_CONSERV_YAC )
     {
@@ -869,7 +863,6 @@ void remap_grids_init(int map_type, bool lextrapolate, int gridID1, remapgrid_t
 
   //if ( src_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
   remap_define_grid(map_type, gridID1, src_grid, "Source");
-
   remap_define_grid(map_type, gridID2, tgt_grid, "Target");
 
   if ( src_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D && tgt_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )

src/util.h

@@ -64,7 +64,7 @@ extern int CDO_opterr;
 extern int CDO_flt_digits;
 extern int CDO_dbl_digits;
 
-extern int remap_genweights;
+extern bool REMAP_genweights;
 
 extern const char *cdoExpName;
 extern int ompNumThreads;