Commit 08028cd1 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added remap_bilinear_scrip.c

parent b7a9fd45
......@@ -590,6 +590,7 @@ src/readline.c -text
src/realtime.c -text
src/remap.h -text
src/remap_bicubic_scrip.c -text
src/remap_bilinear_scrip.c -text
src/remap_conserv.c -text
src/remap_conserv_scrip.c -text
src/remap_distwgt_scrip.c -text
......
......@@ -277,6 +277,7 @@ cdo_SOURCES += Adisit.c \
remap_conserv_scrip.c \
remap_distwgt_scrip.c \
remap_bicubic_scrip.c \
remap_bilinear_scrip.c \
stdnametable.c \
stdnametable.h \
specspace.c \
......
......@@ -165,17 +165,17 @@ am__cdo_SOURCES_DIST = cdo.c Adisit.c Arith.c Arithc.c Arithdays.c \
remapsort.c remap_scrip_io.c remap_search_latbins.c \
remap_store_link.c remap_store_link.h remap_conserv.c \
remap_conserv_scrip.c remap_distwgt_scrip.c \
remap_bicubic_scrip.c stdnametable.c stdnametable.h \
specspace.c specspace.h statistic.c statistic.h table.c text.c \
text.h timebase.h timer.c userlog.c util.c util.h vinterp.c \
vinterp.h zaxis.c clipping/clipping.c clipping/clipping.h \
clipping/area.c clipping/area.h clipping/ensure_array_size.c \
clipping/ensure_array_size.h clipping/geometry_tools.c \
clipping/geometry.h clipping/grid.h clipping/points.h \
clipping/dep_list.h clipping/grid_cell.c clipping/grid_cell.h \
clipping/intersection.c clipping/utils.c clipping/utils.h \
Magplot.c Magvector.c Maggraph.c template_parser.h \
template_parser.c results_template_parser.h \
remap_bicubic_scrip.c remap_bilinear_scrip.c stdnametable.c \
stdnametable.h specspace.c specspace.h statistic.c statistic.h \
table.c text.c text.h timebase.h timer.c userlog.c util.c \
util.h vinterp.c vinterp.h zaxis.c clipping/clipping.c \
clipping/clipping.h clipping/area.c clipping/area.h \
clipping/ensure_array_size.c clipping/ensure_array_size.h \
clipping/geometry_tools.c clipping/geometry.h clipping/grid.h \
clipping/points.h clipping/dep_list.h clipping/grid_cell.c \
clipping/grid_cell.h clipping/intersection.c clipping/utils.c \
clipping/utils.h Magplot.c Magvector.c Maggraph.c \
template_parser.h template_parser.c results_template_parser.h \
results_template_parser.c magics_template_parser.h \
magics_template_parser.c StringUtilities.h StringUtilities.c \
CdoMagicsMapper.h CdoMagicsMapper.c
......@@ -298,7 +298,8 @@ am_cdo_OBJECTS = cdo-cdo.$(OBJEXT) cdo-Adisit.$(OBJEXT) \
cdo-remap_store_link.$(OBJEXT) cdo-remap_conserv.$(OBJEXT) \
cdo-remap_conserv_scrip.$(OBJEXT) \
cdo-remap_distwgt_scrip.$(OBJEXT) \
cdo-remap_bicubic_scrip.$(OBJEXT) cdo-stdnametable.$(OBJEXT) \
cdo-remap_bicubic_scrip.$(OBJEXT) \
cdo-remap_bilinear_scrip.$(OBJEXT) cdo-stdnametable.$(OBJEXT) \
cdo-specspace.$(OBJEXT) cdo-statistic.$(OBJEXT) \
cdo-table.$(OBJEXT) cdo-text.$(OBJEXT) cdo-timer.$(OBJEXT) \
cdo-userlog.$(OBJEXT) cdo-util.$(OBJEXT) cdo-vinterp.$(OBJEXT) \
......@@ -602,16 +603,17 @@ cdo_SOURCES = cdo.c Adisit.c Arith.c Arithc.c Arithdays.c Arithlat.c \
readline.c realtime.c remap.h remaplib.c remapsort.c \
remap_scrip_io.c remap_search_latbins.c remap_store_link.c \
remap_store_link.h remap_conserv.c remap_conserv_scrip.c \
remap_distwgt_scrip.c remap_bicubic_scrip.c stdnametable.c \
stdnametable.h specspace.c specspace.h statistic.c statistic.h \
table.c text.c text.h timebase.h timer.c userlog.c util.c \
util.h vinterp.c vinterp.h zaxis.c clipping/clipping.c \
clipping/clipping.h clipping/area.c clipping/area.h \
clipping/ensure_array_size.c clipping/ensure_array_size.h \
clipping/geometry_tools.c clipping/geometry.h clipping/grid.h \
clipping/points.h clipping/dep_list.h clipping/grid_cell.c \
clipping/grid_cell.h clipping/intersection.c clipping/utils.c \
clipping/utils.h $(am__append_1)
remap_distwgt_scrip.c remap_bicubic_scrip.c \
remap_bilinear_scrip.c stdnametable.c stdnametable.h \
specspace.c specspace.h statistic.c statistic.h table.c text.c \
text.h timebase.h timer.c userlog.c util.c util.h vinterp.c \
vinterp.h zaxis.c clipping/clipping.c clipping/clipping.h \
clipping/area.c clipping/area.h clipping/ensure_array_size.c \
clipping/ensure_array_size.h clipping/geometry_tools.c \
clipping/geometry.h clipping/grid.h clipping/points.h \
clipping/dep_list.h clipping/grid_cell.c clipping/grid_cell.h \
clipping/intersection.c clipping/utils.c clipping/utils.h \
$(am__append_1)
cdo_CPPFLAGS = -I$(top_srcdir)/libcdi/src
cdo_LDADD = $(top_builddir)/libcdi/src/libcdi.la
cdo_LDFLAGS = $(am__append_2)
......@@ -991,6 +993,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-readline.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-realtime.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remap_bicubic_scrip.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remap_bilinear_scrip.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remap_conserv.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remap_conserv_scrip.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remap_distwgt_scrip.Po@am__quote@
......@@ -4339,6 +4342,20 @@ cdo-remap_bicubic_scrip.obj: remap_bicubic_scrip.c
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o cdo-remap_bicubic_scrip.obj `if test -f 'remap_bicubic_scrip.c'; then $(CYGPATH_W) 'remap_bicubic_scrip.c'; else $(CYGPATH_W) '$(srcdir)/remap_bicubic_scrip.c'; fi`
cdo-remap_bilinear_scrip.o: remap_bilinear_scrip.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT cdo-remap_bilinear_scrip.o -MD -MP -MF $(DEPDIR)/cdo-remap_bilinear_scrip.Tpo -c -o cdo-remap_bilinear_scrip.o `test -f 'remap_bilinear_scrip.c' || echo '$(srcdir)/'`remap_bilinear_scrip.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/cdo-remap_bilinear_scrip.Tpo $(DEPDIR)/cdo-remap_bilinear_scrip.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ $(AM_V_CC)source='remap_bilinear_scrip.c' object='cdo-remap_bilinear_scrip.o' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o cdo-remap_bilinear_scrip.o `test -f 'remap_bilinear_scrip.c' || echo '$(srcdir)/'`remap_bilinear_scrip.c
cdo-remap_bilinear_scrip.obj: remap_bilinear_scrip.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT cdo-remap_bilinear_scrip.obj -MD -MP -MF $(DEPDIR)/cdo-remap_bilinear_scrip.Tpo -c -o cdo-remap_bilinear_scrip.obj `if test -f 'remap_bilinear_scrip.c'; then $(CYGPATH_W) 'remap_bilinear_scrip.c'; else $(CYGPATH_W) '$(srcdir)/remap_bilinear_scrip.c'; fi`
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/cdo-remap_bilinear_scrip.Tpo $(DEPDIR)/cdo-remap_bilinear_scrip.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ $(AM_V_CC)source='remap_bilinear_scrip.c' object='cdo-remap_bilinear_scrip.obj' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o cdo-remap_bilinear_scrip.obj `if test -f 'remap_bilinear_scrip.c'; then $(CYGPATH_W) 'remap_bilinear_scrip.c'; else $(CYGPATH_W) '$(srcdir)/remap_bilinear_scrip.c'; fi`
cdo-stdnametable.o: stdnametable.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT cdo-stdnametable.o -MD -MP -MF $(DEPDIR)/cdo-stdnametable.Tpo -c -o cdo-stdnametable.o `test -f 'stdnametable.c' || echo '$(srcdir)/'`stdnametable.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/cdo-stdnametable.Tpo $(DEPDIR)/cdo-stdnametable.Po
......
......@@ -926,7 +926,7 @@ void *Remap(void *argument)
print_remap_info(operfunc, map_type, &remaps[r].src_grid, &remaps[r].tgt_grid, nmiss1);
if ( map_type == MAP_TYPE_CONSERV ) remap_conserv(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_BILINEAR ) remap_bilin(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_BILINEAR ) remap_bilinear(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_BICUBIC ) remap_bicubic(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_DISTWGT ) remap_distwgt(num_neighbors, &remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_CONSPHERE ) remap_consphere(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
......
......@@ -183,7 +183,7 @@ void remap_laf(double *restrict dst_array, double missval, long dst_size, long n
void remap_sum(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict map_wts,
long num_wts, const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array);
void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_consphere(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
......
......@@ -235,4 +235,4 @@ void remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
}
} /* grid_loop1 */
} /* remap_bicub */
} /* remap_bicubic */
#include "cdo.h"
#include "cdo_int.h"
#include "grid.h"
#include "remap.h"
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* */
/* BILINEAR INTERPOLATION */
/* */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
long find_element(double x, long nelem, const double *restrict array);
int rect_grid_search(long *ii, long *jj, double x, double y, long nxm, long nym, const double *restrict xm, const double *restrict ym);
int grid_search_reg2d_nn(long nx, long ny, int *restrict nbr_add, double *restrict nbr_dist, double plat, double plon,
const double *restrict src_center_lat, const double *restrict src_center_lon)
{
int search_result = 0;
long n, srch_add;
long i;
long ii, jj;
long 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;
long jjf = 0, jjl = ny-1;
if ( plon >= src_center_lon[0] && plon <= src_center_lon[nx-1] )
{
if ( src_center_lat[0] < src_center_lat[ny-1] )
{
if ( plat <= src_center_lat[0] )
{ jjf = 0; jjl = 1; }
else
{ jjf = ny-2; jjl = ny-1; }
}
else
{
if ( plat >= src_center_lat[0] )
{ jjf = 0; jjl = 1; }
else
{ jjf = ny-2; jjl = ny-1; }
}
}
sincoslon = malloc(nx*sizeof(double));
for ( 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 )
{
coslat = coslat_dst*cos(src_center_lat[jj]);
sinlat = sinlat_dst*sin(src_center_lat[jj]);
jjskip = jj > 1 && jj < (ny-2);
for ( ii = 0; ii < nx; ++ii )
{
if ( jjskip && ii > 1 && ii < (nx-2) ) continue;
srch_add = jj*nx + ii;
distance = acos(coslat*sincoslon[ii] + sinlat);
if ( distance < dist_min )
{
for ( n = 0; n < 4; ++n )
{
if ( distance < nbr_dist[n] )
{
for ( i = 3; i > n; --i )
{
nbr_add [i] = nbr_add [i-1];
nbr_dist[i] = nbr_dist[i-1];
}
search_result = -1;
nbr_add [n] = srch_add;
nbr_dist[n] = distance;
dist_min = nbr_dist[3];
break;
}
}
}
}
}
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;
return (search_result);
}
int grid_search_reg2d(remapgrid_t *src_grid, int *restrict src_add, double *restrict src_lats,
double *restrict src_lons, double plat, double plon, const int *restrict src_grid_dims,
const double *restrict src_center_lat, const 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
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;
int lfound;
long n;
long nx, nxm, ny;
long ii, iix, jj;
for ( n = 0; n < 4; ++n ) src_add[n] = 0;
nx = src_grid_dims[0];
ny = src_grid_dims[1];
nxm = nx;
if ( src_grid->is_cyclic ) nxm++;
if ( /*plon < 0 &&*/ plon < src_center_lon[0] ) plon += PI2;
if ( /*plon > PI2 &&*/ plon > src_center_lon[nxm-1] ) plon -= PI2;
lfound = rect_grid_search(&ii, &jj, plon, plat, nxm, ny, src_center_lon, src_center_lat);
if ( lfound )
{
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);
src_add[2] = (jj)*nx+(iix);
src_add[3] = (jj)*nx+(ii-1);
src_lons[0] = src_center_lon[ii-1];
src_lons[1] = src_center_lon[iix];
/* For consistency, we must make sure all lons are in same 2pi interval */
if ( src_lons[0] > PI2 ) src_lons[0] -= PI2;
if ( src_lons[0] < 0 ) src_lons[0] += PI2;
if ( src_lons[1] > PI2 ) src_lons[1] -= PI2;
if ( src_lons[1] < 0 ) src_lons[1] += PI2;
src_lons[2] = src_lons[1];
src_lons[3] = src_lons[0];
src_lats[0] = src_center_lat[jj-1];
src_lats[1] = src_lats[0];
src_lats[2] = src_center_lat[jj];
src_lats[3] = src_lats[2];
search_result = 1;
return (search_result);
}
/*
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 average for this point\n");
*/
search_result = grid_search_reg2d_nn(nx, ny, src_add, src_lats, plat, plon, src_center_lat, src_center_lon);
return (search_result);
} /* grid_search_reg2d */
static
int grid_search_nn(long min_add, long max_add, int *restrict nbr_add, double *restrict nbr_dist,
double plat, double plon,
const double *restrict src_center_lat, const double *restrict src_center_lon)
{
int search_result = 0;
long n, srch_add;
long i;
double dist_min, distance; /* For computing dist-weighted avg */
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;
for ( srch_add = min_add; srch_add <= max_add; ++srch_add )
{
distance = acos(coslat_dst*cos(src_center_lat[srch_add])*
(coslon_dst*cos(src_center_lon[srch_add]) +
sinlon_dst*sin(src_center_lon[srch_add]))+
sinlat_dst*sin(src_center_lat[srch_add]));
if ( distance < dist_min )
{
for ( n = 0; n < 4; ++n )
{
if ( distance < nbr_dist[n] )
{
for ( i = 3; i > n; --i )
{
nbr_add [i] = nbr_add [i-1];
nbr_dist[i] = nbr_dist[i-1];
}
search_result = -1;
nbr_add [n] = srch_add;
nbr_dist[n] = distance;
dist_min = nbr_dist[3];
break;
}
}
}
}
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;
return (search_result);
}
int grid_search(remapgrid_t *src_grid, int *restrict src_add, double *restrict src_lats,
double *restrict src_lons, double plat, double plon, const int *restrict src_grid_dims,
const double *restrict src_center_lat, const double *restrict src_center_lon,
const restr_t *restrict src_grid_bound_box, const int *restrict src_bin_add)
{
/*
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
double src_center_lat[] ! latitude of each src grid center
double src_center_lon[] ! longitude of each src grid center
restr_t src_grid_bound_box[][4] ! bound box for source grid
int src_bin_add[][2] ! latitude bins for restricting
*/
/* Local variables */
long n, n2, next_n, srch_add, srch_add4; /* dummy indices */
long nx, ny; /* dimensions of src grid */
long min_add, max_add; /* addresses for restricting search */
long i, j, jp1, ip1, n_add, e_add, ne_add; /* addresses */
long nbins;
/* Vectors for cross-product check */
double vec1_lat, vec1_lon;
double vec2_lat, vec2_lon, cross_product;
int scross[4], scross_last = 0;
int search_result = 0;
restr_t rlat, rlon;
restr_t *bin_lats = src_grid->bin_lats;
nbins = src_grid->num_srch_bins;
rlat = RESTR_SCALE(plat);
rlon = RESTR_SCALE(plon);
/* restrict search first using bins */
for ( n = 0; n < 4; ++n ) src_add[n] = 0;
min_add = src_grid->size-1;
max_add = 0;
for ( n = 0; n < nbins; ++n )
{
n2 = n<<1;
if ( rlat >= bin_lats[n2] && rlat <= bin_lats[n2+1] )
{
if ( src_bin_add[n2 ] < min_add ) min_add = src_bin_add[n2 ];
if ( src_bin_add[n2+1] > max_add ) max_add = src_bin_add[n2+1];
}
}
/* Now perform a more detailed search */
nx = src_grid_dims[0];
ny = src_grid_dims[1];
/* srch_loop */
for ( srch_add = min_add; srch_add <= max_add; ++srch_add )
{
srch_add4 = srch_add<<2;
/* First check bounding box */
if ( rlon >= src_grid_bound_box[srch_add4+2] &&
rlon <= src_grid_bound_box[srch_add4+3] &&
rlat >= src_grid_bound_box[srch_add4 ] &&
rlat <= src_grid_bound_box[srch_add4+1])
{
/* We are within bounding box so get really serious */
/* Determine neighbor addresses */
j = srch_add/nx;
i = srch_add - j*nx;
if ( i < (nx-1) )
ip1 = i + 1;
else
{
/* 2009-01-09 Uwe Schulzweida: bug fix */
if ( src_grid->is_cyclic )
ip1 = 0;
else
ip1 = i;
}
if ( j < (ny-1) )
jp1 = j + 1;
else
{
/* 2008-12-17 Uwe Schulzweida: latitute cyclic ??? (bug fix) */
jp1 = j;
}
n_add = jp1*nx + i;
e_add = j *nx + ip1;
ne_add = jp1*nx + ip1;
src_lons[0] = src_center_lon[srch_add];
src_lons[1] = src_center_lon[e_add];
src_lons[2] = src_center_lon[ne_add];
src_lons[3] = src_center_lon[n_add];
src_lats[0] = src_center_lat[srch_add];
src_lats[1] = src_center_lat[e_add];
src_lats[2] = src_center_lat[ne_add];
src_lats[3] = src_center_lat[n_add];
/* For consistency, we must make sure all lons are in same 2pi interval */
vec1_lon = src_lons[0] - plon;
if ( vec1_lon > PI ) src_lons[0] -= PI2;
else if ( vec1_lon < -PI ) src_lons[0] += PI2;
for ( n = 1; n < 4; ++n )
{
vec1_lon = src_lons[n] - src_lons[0];
if ( vec1_lon > PI ) src_lons[n] -= PI2;
else if ( vec1_lon < -PI ) src_lons[n] += PI2;
}
/* corner_loop */
for ( n = 0; n < 4; ++n )
{
next_n = (n+1)%4;
/*
Here we take the cross product of the vector making
up each box side with the vector formed by the vertex
and search point. If all the cross products are
positive, the point is contained in the box.
*/
vec1_lat = src_lats[next_n] - src_lats[n];
vec1_lon = src_lons[next_n] - src_lons[n];
vec2_lat = plat - src_lats[n];
vec2_lon = plon - src_lons[n];
/* Check for 0,2pi crossings */
if ( vec1_lon > THREE*PIH ) vec1_lon -= PI2;
else if ( vec1_lon < -THREE*PIH ) vec1_lon += PI2;
if ( vec2_lon > THREE*PIH ) vec2_lon -= PI2;
else if ( vec2_lon < -THREE*PIH ) vec2_lon += PI2;
cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat;
/* If cross product is less than ZERO, this cell doesn't work */
/* 2008-10-16 Uwe Schulzweida: bug fix for cross_product eq zero */
scross[n] = cross_product < 0 ? -1 : cross_product > 0 ? 1 : 0;
if ( n == 0 ) scross_last = scross[n];
if ( (scross[n] < 0 && scross_last > 0) || (scross[n] > 0 && scross_last < 0) ) break;
scross_last = scross[n];
} /* corner_loop */
if ( n >= 4 )
{
n = 0;
if ( scross[0]>=0 && scross[1]>=0 && scross[2]>=0 && scross[3]>=0 ) n = 4;
else if ( scross[0]<=0 && scross[1]<=0 && scross[2]<=0 && scross[3]<=0 ) n = 4;
}
/* If cross products all same sign, we found the location */
if ( n >= 4 )
{
src_add[0] = srch_add;
src_add[1] = e_add;
src_add[2] = ne_add;
src_add[3] = n_add;
search_result = 1;
return (search_result);
}
/* Otherwise move on to next cell */
} /* Bounding box check */
} /* srch_loop */
/*
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 average for this point\n");
*/
search_result = grid_search_nn(min_add, max_add, src_add, src_lats, plat, plon, src_center_lat, src_center_lon);
return (search_result);
} /* grid_search */
/*
This routine stores the address and weight for four links associated with one destination
point in the appropriate address and weight arrays and resizes those arrays if necessary.
*/
void store_link_bilin(remapvars_t *rv, int dst_add, int *restrict src_add, double *restrict weights)
{
/*
Input variables:
int dst_add ! address on destination grid
int src_add[4] ! addresses on source grid
double weights[4] ! array of remapping weights for these links
*/
long n;
long nlink;
/*
Increment number of links and check to see if remap arrays need
to be increased to accomodate the new link. Then store the link.