added remapgrid_get_lonlat

src/remap.h

@@ -234,5 +234,6 @@ int grid_search(remapgrid_t *src_grid, int *restrict src_add, double *restrict s
 int find_ij_weights(double plon, double plat, double* restrict src_lats, double* restrict src_lons, double *ig, double *jg);
 int rect_grid_search(long *ii, long *jj, double x, double y, long nxm, long nym, const double *restrict xm, const double *restrict ym);
 
+void remapgrid_get_lonlat(remapgrid_t *grid, unsigned cell_add, double *plon, double *plat);
 
 #endif  /* _REMAP_H */

src/remap_distwgt.c

@@ -448,7 +448,7 @@ void distwgt_remap(double* restrict tgt_point, const double* restrict src_array,
 void remap_distwgt(unsigned num_neighbors, remapgrid_t *src_grid, remapgrid_t *tgt_grid, const double* restrict src_array, double* restrict tgt_array, double missval)
 {
   /*  Local variables */
-  int remap_grid_type = src_grid->remap_grid_type;
+  int src_remap_grid_type = src_grid->remap_grid_type;
 
   if ( cdoVerbose ) cdoPrint("Called %s()", __func__);
 
@@ -470,7 +470,7 @@ void remap_distwgt(unsigned num_neighbors, remapgrid_t *src_grid, remapgrid_t *t
   start = clock();
 
   struct gridsearch *gs = NULL;
-  if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
+  if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
     gs = gridsearch_create_reg2d(nx, ny, 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);
@@ -487,7 +487,7 @@ void remap_distwgt(unsigned num_neighbors, remapgrid_t *src_grid, remapgrid_t *t
 
 #if defined(_OPENMP)
 #pragma omp parallel for default(none) \
-  shared(gs, num_neighbors, remap_grid_type, src_grid, tgt_grid, tgt_grid_size, findex) \
+  shared(gs, num_neighbors, src_remap_grid_type, src_grid, tgt_grid, tgt_grid_size, findex) \
   shared(src_array, tgt_array, missval) \
   private(nbr_mask, nbr_add, nbr_dist)
 #endif
@@ -502,12 +502,12 @@ void remap_distwgt(unsigned num_neighbors, remapgrid_t *src_grid, remapgrid_t *t
       tgt_array[tgt_cell_add] = missval;
 
       if ( ! tgt_grid->mask[tgt_cell_add] ) continue;
-	
-      double plat = tgt_grid->cell_center_lat[tgt_cell_add];
-      double plon = tgt_grid->cell_center_lon[tgt_cell_add];
+
+      double plon = 0, plat = 0;
+      remapgrid_get_lonlat(tgt_grid, tgt_cell_add, &plon, &plat);
 
       /* Find nearest grid points on source grid and distances to each point */
-      if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
+      if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
 	grid_search_nbr_reg2d(gs, num_neighbors, src_grid, nbr_add, nbr_dist, 
 			      plat, plon, src_grid->dims);
       else

src/remaplib.c

@@ -371,6 +371,25 @@ void check_lon_range2(long nc, long nlons, double *corners, double *centers)
     }
 }
 
+
+void remapgrid_get_lonlat(remapgrid_t *grid, unsigned cell_add, double *plon, double *plat)
+{
+  if ( grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )
+    {
+      unsigned nx = grid->dims[0];
+      unsigned iy = cell_add/nx;
+      unsigned ix = cell_add - iy*nx;
+      *plat = grid->reg2d_center_lat[iy];
+      *plon = grid->reg2d_center_lon[ix];
+      if ( *plon < 0 ) *plon += PI2;
+    }
+  else
+    {
+      *plat = grid->cell_center_lat[cell_add];
+      *plon = grid->cell_center_lon[cell_add];
+    }
+}
+
 static
 void check_lon_range(long nlons, double *lons)
 {
@@ -647,6 +666,7 @@ void remap_define_grid(int map_type, int gridID, remapgrid_t *grid)
   int lgrid_gen_bounds = FALSE;
   int gridID_gme = -1;
 
+  printf("gridID %d %d %d\n", gridID, remap_write_remap == TRUE, grid->remap_grid_type != REMAP_GRID_TYPE_REG2D);
   if ( gridInqType(grid->gridID) != GRID_UNSTRUCTURED && gridInqType(grid->gridID) != GRID_CURVILINEAR )
     {
       if ( gridInqType(grid->gridID) == GRID_GME )
@@ -660,7 +680,9 @@ void remap_define_grid(int map_type, int gridID, remapgrid_t *grid)
       else if ( remap_write_remap == TRUE || grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
 	{
 	  lgrid_destroy = TRUE;
+          printf("to curvilinear %d\n", gridID);
 	  gridID = gridToCurvilinear(grid->gridID, 1);
+          printf("   curvilinear %d\n", gridID);
 	  lgrid_gen_bounds = TRUE;
 	}
     }
@@ -830,6 +852,7 @@ void cell_bounding_boxes(remapgrid_t *grid, int remap_grid_basis)
 
 void remap_grids_init(int map_type, int lextrapolate, int gridID1, remapgrid_t *src_grid, int gridID2, remapgrid_t *tgt_grid)
 {
+  int lbounds = TRUE;
   int reg2d_src_gridID = gridID1;
   int reg2d_tgt_gridID = gridID2;
 
@@ -844,16 +867,15 @@ void remap_grids_init(int map_type, int lextrapolate, int gridID1, remapgrid_t *
       // src_grid->remap_grid_type = 0;
     }
 
-  if ( src_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D && map_type == MAP_TYPE_CONSERV_YAC )
+  if ( src_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )
     {
-      if ( IS_REG2D_GRID(gridID2) ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
+      if ( IS_REG2D_GRID(gridID2) && map_type == MAP_TYPE_CONSERV_YAC ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
       // else src_grid->remap_grid_type = -1;
     }
 
-  if ( remap_gen_weights == FALSE && map_type == MAP_TYPE_BILINEAR )
+  if ( remap_gen_weights == FALSE && IS_REG2D_GRID(gridID2) && tgt_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
     {
-      if ( IS_REG2D_GRID(gridID2) ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
-      // else src_grid->remap_grid_type = -1;
+      if ( map_type == MAP_TYPE_BILINEAR || map_type == MAP_TYPE_DISTWGT ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
     }
 
   if ( lextrapolate > 0 )
@@ -917,7 +939,7 @@ void remap_grids_init(int map_type, int lextrapolate, int gridID1, remapgrid_t *
 	gridInqType(src_grid->gridID) == GRID_LAEA || 
 	gridInqType(src_grid->gridID) == GRID_SINUSOIDAL) )
     {
-      src_grid->gridID = gridID1 = gridToCurvilinear(src_grid->gridID, 1);
+      src_grid->gridID = gridID1 = gridToCurvilinear(src_grid->gridID, lbounds);
     }
 
   if ( !src_grid->lextrapolate && gridInqSize(src_grid->gridID) > 1 &&