Commit 8843484b authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

TEST_KDTREE

parent f7496e69
......@@ -202,6 +202,83 @@ void scrip_remap_bicubic_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
-----------------------------------------------------------------------
*/
//#define TEST_KDTREE
#ifdef TEST_KDTREE
#include "grid_search.h"
int grid_search_test(struct gridsearch *gs, size_t *restrict src_add, double *restrict src_lats,
double *restrict src_lons, double plat, double plon, const size_t *restrict src_grid_dims,
double *restrict src_center_lat, double *restrict src_center_lon)
{
/*
Output variables:
int src_add[4] ! address of each corner point enclosing P
double src_lats[4] ! latitudes of the four corner points
double src_lons[4] ! longitudes of the four corner points
Input variables:
double plat ! latitude of the search point
double plon ! longitude of the search point
int src_grid_dims[2] ! size of each src grid dimension
*/
bool is_cyclic = true;
int search_result = 0;
for ( unsigned n = 0; n < 4; ++n ) src_add[n] = 0;
/* Now perform a more detailed search */
size_t nx = src_grid_dims[0];
size_t ny = src_grid_dims[1];
double search_radius = gs->search_radius;
const double range0 = SQR(search_radius);
double range = range0;
size_t add = gridsearch_nearest(gs, plon, plat, &range);
// printf("plon, plat, add, range %g %g %g %g %zu %g\n", plon*RAD2DEG, plat*RAD2DEG,
// src_center_lon[add]*RAD2DEG, src_center_lat[add]*RAD2DEG,add, range);
if ( add != GS_NOT_FOUND )
{
size_t idx[4];
for ( unsigned k = 0; k < 4; ++k )
{
/* Determine neighbor addresses */
size_t j = add/nx;
size_t i = add - 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, src_add, src_lons, src_lats,
plon, plat, src_center_lon, src_center_lat) )
{
search_result = 1;
return search_result;
}
/* Otherwise move on to next cell */
}
/*
If no cell found, point is likely either in a box that straddles either pole or is outside
the grid. Fall back to a distance-weighted average of the four closest points.
Go ahead and compute weights here, but store in src_lats and return -add to prevent the
parent routine from computing bilinear weights.
*/
// if ( !src_grid->lextrapolate ) return search_result;
/*
printf("Could not find location for %g %g\n", plat*RAD2DEG, plon*RAD2DEG);
printf("Using nearest-neighbor for this point\n");
*/
search_result = add;
}
return search_result;
} /* grid_search_test */
#endif
void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const double* restrict src_array, double* restrict tgt_array, double missval)
{
int remap_grid_type = src_grid->remap_grid_type;
......@@ -210,6 +287,14 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
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(xIsCyclic, dims, src_grid->size, src_grid->cell_center_lon, src_grid->cell_center_lat);
#endif
size_t tgt_grid_size = tgt_grid->size;
/* Compute mappings from source to target grid */
......@@ -258,10 +343,15 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
plat, plon, src_grid->dims,
src_grid->reg2d_center_lat, src_grid->reg2d_center_lon);
else
#ifdef TEST_KDTREE
search_result = grid_search_test(gs, src_add, src_lats, src_lons, plat, plon, src_grid->dims,
src_grid->cell_center_lat, src_grid->cell_center_lon);
#else
search_result = grid_search(src_grid, src_add, src_lats, src_lons,
plat, plon, src_grid->dims,
src_grid->cell_center_lat, src_grid->cell_center_lon,
src_grid->cell_bound_box, src_grid->bin_addr);
#endif
/* Check to see if points are land points */
if ( search_result > 0 )
......@@ -316,4 +406,8 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
Free(grad1_lon);
Free(grad1_latlon);
} /* scrip_remap_bicubic */
#ifdef TEST_KDTREE
if ( gs ) gridsearch_delete(gs);
#endif
} // scrip_remap_bicubic
......@@ -374,7 +374,7 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
bool xIsCyclic = false;
size_t dims[2] = {src_grid->size, 0};
struct gridsearch *gs = NULL;
if (remap_grid_type != REMAP_GRID_TYPE_REG2D )
if ( remap_grid_type != REMAP_GRID_TYPE_REG2D )
gs = gridsearch_create_nn(xIsCyclic, dims, src_grid->size, src_grid->cell_center_lon, src_grid->cell_center_lat);
#endif
......@@ -391,7 +391,7 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
#if defined(_OPENMP)
#pragma omp parallel for default(none) schedule(static) \
shared(cdoSilentMode, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, findex)
shared(cdoSilentMode, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, findex)
#endif
for ( size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
{
......
......@@ -799,7 +799,9 @@ void remap_grids_init(int map_type, bool lextrapolate, int gridID1, remapgrid_t
{
remap_define_reg2d(reg2d_tgt_gridID, tgt_grid);
}
else if ( map_type != MAP_TYPE_DISTWGT )
else if ( map_type != MAP_TYPE_DISTWGT
// && map_type != MAP_TYPE_BILINEAR
)
{
cell_bounding_boxes(src_grid, REMAP_GRID_BASIS_SRC);
cell_bounding_boxes(tgt_grid, REMAP_GRID_BASIS_TGT);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment