Skip to content
GitLab
Commit b7a9fd45 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added remap_bicubic_scrip.c

parent a60df6cc
Hide whitespace changes
Inline Side-by-side
.gitattributes
View file @ b7a9fd45
......@@ -589,6 +589,7 @@ src/pthread_debug.h -text
src/readline.c -text
src/realtime.c -text
src/remap.h -text
src/remap_bicubic_scrip.c -text
src/remap_conserv.c -text
src/remap_conserv_scrip.c -text
src/remap_distwgt_scrip.c -text
......
src/Makefile.am
View file @ b7a9fd45
......@@ -276,6 +276,7 @@ cdo_SOURCES += Adisit.c \
remap_conserv.c \
remap_conserv_scrip.c \
remap_distwgt_scrip.c \
remap_bicubic_scrip.c \
stdnametable.c \
stdnametable.h \
specspace.c \
......
src/Makefile.in
View file @ b7a9fd45
......@@ -164,17 +164,18 @@ am__cdo_SOURCES_DIST = cdo.c Adisit.c Arith.c Arithc.c Arithdays.c \
pthread_debug.h 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 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_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 \
results_template_parser.c magics_template_parser.h \
magics_template_parser.c StringUtilities.h StringUtilities.c \
CdoMagicsMapper.h CdoMagicsMapper.c
......@@ -296,7 +297,8 @@ am_cdo_OBJECTS = cdo-cdo.$(OBJEXT) cdo-Adisit.$(OBJEXT) \
cdo-remap_search_latbins.$(OBJEXT) \
cdo-remap_store_link.$(OBJEXT) cdo-remap_conserv.$(OBJEXT) \
cdo-remap_conserv_scrip.$(OBJEXT) \
cdo-remap_distwgt_scrip.$(OBJEXT) cdo-stdnametable.$(OBJEXT) \
cdo-remap_distwgt_scrip.$(OBJEXT) \
cdo-remap_bicubic_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) \
......@@ -600,16 +602,16 @@ 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 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 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)
......@@ -988,6 +990,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-pthread_debug.Po@am__quote@
@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_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@
......@@ -4322,6 +4325,20 @@ cdo-remap_distwgt_scrip.obj: remap_distwgt_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_distwgt_scrip.obj `if test -f 'remap_distwgt_scrip.c'; then $(CYGPATH_W) 'remap_distwgt_scrip.c'; else $(CYGPATH_W) '$(srcdir)/remap_distwgt_scrip.c'; fi`
cdo-remap_bicubic_scrip.o: remap_bicubic_scrip.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT cdo-remap_bicubic_scrip.o -MD -MP -MF $(DEPDIR)/cdo-remap_bicubic_scrip.Tpo -c -o cdo-remap_bicubic_scrip.o `test -f 'remap_bicubic_scrip.c' || echo '$(srcdir)/'`remap_bicubic_scrip.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/cdo-remap_bicubic_scrip.Tpo $(DEPDIR)/cdo-remap_bicubic_scrip.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ $(AM_V_CC)source='remap_bicubic_scrip.c' object='cdo-remap_bicubic_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_bicubic_scrip.o `test -f 'remap_bicubic_scrip.c' || echo '$(srcdir)/'`remap_bicubic_scrip.c
cdo-remap_bicubic_scrip.obj: remap_bicubic_scrip.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT cdo-remap_bicubic_scrip.obj -MD -MP -MF $(DEPDIR)/cdo-remap_bicubic_scrip.Tpo -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`
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/cdo-remap_bicubic_scrip.Tpo $(DEPDIR)/cdo-remap_bicubic_scrip.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ $(AM_V_CC)source='remap_bicubic_scrip.c' object='cdo-remap_bicubic_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_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-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
......
src/Remap.c
View file @ b7a9fd45
......@@ -927,7 +927,7 @@ void *Remap(void *argument)
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_BICUBIC ) remap_bicub(&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);
......
src/YAR.c
View file @ b7a9fd45
......@@ -139,77 +139,6 @@ void set_source_data(double * source_data, double init_value,
}
static
long find_ij_weights(double plon, double plat, double *restrict src_lats, double *restrict src_lons, double *ig, double *jg)
{
#define THREE 3.0
#define HALF 0.5
long Max_Iter = 100;
double converge = 1.e-10; /* Convergence criterion */
long iter; /* iteration counters */
double iguess, jguess; /* current guess for bilinear coordinate */
double deli, delj; /* corrections to i,j */
double dth1, dth2, dth3; /* some latitude differences */
double dph1, dph2, dph3; /* some longitude differences */
double dthp, dphp; /* difference between point and sw corner */
double mat1, mat2, mat3, mat4; /* matrix elements */
double determinant; /* matrix determinant */
/* Iterate to find i,j for bilinear approximation */
dth1 = src_lats[1] - src_lats[0];
dth2 = src_lats[3] - src_lats[0];
dth3 = src_lats[2] - src_lats[1] - dth2;
dph1 = src_lons[1] - src_lons[0];
dph2 = src_lons[3] - src_lons[0];
dph3 = src_lons[2] - src_lons[1];
if ( dph1 > THREE*PIH ) dph1 -= PI2;
if ( dph2 > THREE*PIH ) dph2 -= PI2;
if ( dph3 > THREE*PIH ) dph3 -= PI2;
if ( dph1 < -THREE*PIH ) dph1 += PI2;
if ( dph2 < -THREE*PIH ) dph2 += PI2;
if ( dph3 < -THREE*PIH ) dph3 += PI2;
dph3 = dph3 - dph2;
iguess = HALF;
jguess = HALF;
for ( iter = 0; iter < Max_Iter; ++iter )
{
dthp = plat - src_lats[0] - dth1*iguess - dth2*jguess - dth3*iguess*jguess;
dphp = plon - src_lons[0];
if ( dphp > THREE*PIH ) dphp -= PI2;
if ( dphp < -THREE*PIH ) dphp += PI2;
dphp = dphp - dph1*iguess - dph2*jguess - dph3*iguess*jguess;
mat1 = dth1 + dth3*jguess;
mat2 = dth2 + dth3*iguess;
mat3 = dph1 + dph3*jguess;
mat4 = dph2 + dph3*iguess;
determinant = mat1*mat4 - mat2*mat3;
deli = (dthp*mat4 - dphp*mat2)/determinant;
delj = (dphp*mat1 - dthp*mat3)/determinant;
if ( fabs(deli) < converge && fabs(delj) < converge ) break;
iguess += deli;
jguess += delj;
}
*ig = iguess;
*jg = jguess;
return (iter);
}
void yar_remap_bil(field_t *field1, field_t *field2)
{
int nlonIn, nlatIn;
......
src/remap.h
View file @ b7a9fd45
......@@ -184,7 +184,7 @@ void remap_sum(double *restrict dst_array, double missval, long dst_size, long n
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_bicub(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);
void remap_distwgt(int num_neighbors, remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
......@@ -217,5 +217,17 @@ long get_srch_cells(long tgt_grid_add, long nbins, int *bin_addr1, int *bin_addr
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 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);
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);
long 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);
#endif /* _REMAP_H */
src/remap_bicubic_scrip.c 0 → 100644
View file @ b7a9fd45
#include "cdo.h"
#include "cdo_int.h"
#include "grid.h"
#include "remap.h"
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* */
/* BICUBIC INTERPOLATION */
/* */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/*
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.
*/
static
void store_link_bicub(remapvars_t *rv, int dst_add, int *restrict src_add, double weights[4][4])
{
/*
Input variables:
int dst_add ! address on destination grid
int src_add[4] ! addresses on source grid
double weights[4][4] ! array of remapping weights for these links
*/
long n, k;
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.
*/
nlink = rv->num_links;
rv->num_links += 4;
if ( rv->num_links >= rv->max_links )
resize_remap_vars(rv, rv->resize_increment);
if ( rv->sort_add == FALSE )
{
for ( n = 0; n < 3; ++n )
if ( src_add[n] > src_add[n+1] ) break;
if ( n < 2 ) rv->sort_add = TRUE;
else if ( n == 2 ) // swap 3rd and 4th elements
{
{
int itmp;
itmp = src_add[2];
src_add[2] = src_add[3];
src_add[3] = itmp;
}
{
double dtmp[4];
for ( k = 0; k < 4; ++k ) dtmp[k] = weights[k][2];
for ( k = 0; k < 4; ++k ) weights[k][2] = weights[k][3];
for ( k = 0; k < 4; ++k ) weights[k][3] = dtmp[k];
}
}
}
for ( n = 0; n < 4; ++n )
{
rv->src_grid_add[nlink+n] = src_add[n];
rv->tgt_grid_add[nlink+n] = dst_add;
for ( k = 0; k < 4; ++k )
rv->wts[4*(nlink+n)+k] = weights[k][n];
}
} /* store_link_bicub */
/*
-----------------------------------------------------------------------
This routine computes the weights for a bicubic interpolation.
-----------------------------------------------------------------------
*/
void remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
{
/* Local variables */
int lwarn = TRUE;
int search_result;
long tgt_grid_size;
long n, icount;
long dst_add; /* destination addresss */
long iter; /* iteration counters */
int src_add[4]; /* address for the four source points */
double src_lats[4]; /* latitudes of four bilinear corners */
double src_lons[4]; /* longitudes of four bilinear corners */
double wgts[4][4]; /* bicubic weights for four corners */
double plat, plon; /* lat/lon coords of destination point */
double findex = 0;
extern long remap_max_iter;
int remap_grid_type = src_grid->remap_grid_type;
if ( cdoVerbose ) cdoPrint("Called %s()", __func__);
progressInit();
tgt_grid_size = tgt_grid->size;
/* Compute mappings from source to target grid */
if ( src_grid->rank != 2 )
cdoAbort("Can not do bicubic interpolation when source grid rank != 2");
/* Loop over destination grid */
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, remap_max_iter, lwarn, findex) \
private(dst_add, n, icount, iter, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
schedule(dynamic,1)
#endif
/* grid_loop1 */
for ( dst_add = 0; dst_add < tgt_grid_size; ++dst_add )
{
int lprogress = 1;
#if defined(_OPENMP)
if ( omp_get_thread_num() != 0 ) lprogress = 0;
#endif
#if defined(_OPENMP)
#pragma omp atomic
#endif
findex++;
if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);
if ( ! tgt_grid->mask[dst_add] ) continue;
plat = tgt_grid->cell_center_lat[dst_add];
plon = tgt_grid->cell_center_lon[dst_add];
/* Find nearest square of grid points on source grid */
if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
search_result = grid_search_reg2d(src_grid, src_add, src_lats, src_lons,
plat, plon, src_grid->dims,
src_grid->reg2d_center_lat, src_grid->reg2d_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);
/* Check to see if points are land points */
if ( search_result > 0 )
{
for ( n = 0; n < 4; ++n )
if ( ! src_grid->mask[src_add[n]] ) search_result = 0;
}
/* If point found, find local iw,jw coordinates for weights */
if ( search_result > 0 )
{
double iw, jw; /* current guess for bilinear coordinate */
tgt_grid->cell_frac[dst_add] = ONE;
iter = find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw);
if ( iter < remap_max_iter )
{
/* Successfully found iw,jw - compute weights */
wgts[0][0] = (ONE-jw*jw*(THREE-TWO*jw)) * (ONE-iw*iw*(THREE-TWO*iw));
wgts[0][1] = (ONE-jw*jw*(THREE-TWO*jw)) * iw*iw*(THREE-TWO*iw);
wgts[0][2] = jw*jw*(THREE-TWO*jw) * iw*iw*(THREE-TWO*iw);
wgts[0][3] = jw*jw*(THREE-TWO*jw) * (ONE-iw*iw*(THREE-TWO*iw));
wgts[1][0] = (ONE-jw*jw*(THREE-TWO*jw)) * iw*(iw-ONE)*(iw-ONE);
wgts[1][1] = (ONE-jw*jw*(THREE-TWO*jw)) * iw*iw*(iw-ONE);
wgts[1][2] = jw*jw*(THREE-TWO*jw) * iw*iw*(iw-ONE);
wgts[1][3] = jw*jw*(THREE-TWO*jw) * iw*(iw-ONE)*(iw-ONE);
wgts[2][0] = jw*(jw-ONE)*(jw-ONE) * (ONE-iw*iw*(THREE-TWO*iw));
wgts[2][1] = jw*(jw-ONE)*(jw-ONE) * iw*iw*(THREE-TWO*iw);
wgts[2][2] = jw*jw*(jw-ONE) * iw*iw*(THREE-TWO*iw);
wgts[2][3] = jw*jw*(jw-ONE) * (ONE-iw*iw*(THREE-TWO*iw));
wgts[3][0] = iw*(iw-ONE)*(iw-ONE) * jw*(jw-ONE)*(jw-ONE);
wgts[3][1] = iw*iw*(iw-ONE) * jw*(jw-ONE)*(jw-ONE);
wgts[3][2] = iw*iw*(iw-ONE) * jw*jw*(jw-ONE);
wgts[3][3] = iw*(iw-ONE)*(iw-ONE) * jw*jw*(jw-ONE);
#if defined(_OPENMP)
#pragma omp critical
#endif
store_link_bicub(rv, dst_add, src_add, wgts);
}
else
{
if ( cdoVerbose || lwarn )
{
lwarn = FALSE;
// cdoWarning("Iteration for iw,jw exceed max iteration count of %d!", remap_max_iter);
cdoWarning("Bicubic interpolation failed for some grid points - used a distance-weighted average instead!");
}
search_result = -1;
}
}
/*
Search for bicubic failed - use a distance-weighted average instead (this is typically near the pole)
*/
if ( search_result < 0 )
{
icount = 0;
for ( n = 0; n < 4; ++n )
{
if ( src_grid->mask[src_add[n]] )
icount++;
else
src_lats[n] = ZERO;
}
if ( icount > 0 )
{
/* Renormalize weights */
double sum_wgts = 0.0; /* sum of weights for normalization */
/* 2012-05-08 Uwe Schulzweida: using absolute value of src_lats (bug fix) */
for ( n = 0; n < 4; ++n ) sum_wgts += fabs(src_lats[n]);
for ( n = 0; n < 4; ++n ) wgts[0][n] = fabs(src_lats[n])/sum_wgts;
for ( n = 0; n < 4; ++n ) wgts[1][n] = ZERO;
for ( n = 0; n < 4; ++n ) wgts[2][n] = ZERO;
for ( n = 0; n < 4; ++n ) wgts[3][n] = ZERO;
tgt_grid->cell_frac[dst_add] = ONE;
#if defined(_OPENMP)
#pragma omp critical
#endif
store_link_bicub(rv, dst_add, src_add, wgts);
}
}
} /* grid_loop1 */
} /* remap_bicub */
src/remaplib.c
View file @ b7a9fd45
......@@ -68,7 +68,7 @@
static int remap_write_remap = FALSE;
static int remap_num_srch_bins = 180;
#define DEFAULT_MAX_ITER 100
static long remap_max_iter = DEFAULT_MAX_ITER; /* Max iteration count for i, j iteration */
long remap_max_iter = DEFAULT_MAX_ITER; /* Max iteration count for i, j iteration */
void remap_set_int(int remapvar, int value)
{
......@@ -1428,9 +1428,6 @@ void remap_sum(double *restrict dst_array, double missval, long dst_size, long n
}
static double converge = 1.e-10; /* Convergence criterion */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* */
/* BILINEAR INTERPOLATION */
......@@ -1531,7 +1528,7 @@ int grid_search_reg2d_nn(long nx, long ny, int *restrict nbr_add, double *restri
return (search_result);
}
static
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)
......@@ -1672,7 +1669,7 @@ int grid_search_nn(long min_add, long max_add, int *restrict nbr_add, double *re
return (search_result);
}
static
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,
......@@ -1939,8 +1936,8 @@ void store_link_bilin(remapvars_t *rv, int dst_add, int *restrict src_add, doubl
} /* store_link_bilin */
static
long find_ij_weights(double plon, double plat, double *restrict src_lats, double *restrict src_lons, double *ig, double *jg)
long find_ij_weights(double plon, double plat, double* restrict src_lats, double* restrict src_lons, double *ig, double *jg)
{
long iter; /* iteration counters */
double iguess, jguess; /* current guess for bilinear coordinate */
......@@ -1950,6 +1947,7 @@ long find_ij_weights(double plon, double plat, double *restrict src_lats, double
double dthp, dphp; /* difference between point and sw corner */
double mat1, mat2, mat3, mat4; /* matrix elements */
double determinant; /* matrix determinant */
double converge = 1.e-10; /* Convergence criterion */
/* Iterate to find iw,jw for bilinear approximation */
......@@ -2049,7 +2047,7 @@ void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, remap_max_iter, converge, lwarn, findex) \
shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, remap_max_iter, lwarn, findex) \
private(dst_add, n, icount, iter, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
schedule(dynamic,1)
#endif
......@@ -2177,237 +2175,6 @@ void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
} /* remap_bilin */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* */
/* BICUBIC INTERPOLATION */
/* */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/*
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.
*/
static
void store_link_bicub(remapvars_t *rv, int dst_add, int *restrict src_add, double weights[4][4])
{
/*
Input variables:
int dst_add ! address on destination grid
int src_add[4] ! addresses on source grid
double weights[4][4] ! array of remapping weights for these links
*/
long n, k;
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.
*/
nlink = rv->num_links;
rv->num_links += 4;
if ( rv->num_links >= rv->max_links )
resize_remap_vars(rv, rv->resize_increment);
if ( rv->sort_add == FALSE )
{
for ( n = 0; n < 3; ++n )
if ( src_add[n] > src_add[n+1] ) break;
if ( n < 2 ) rv->sort_add = TRUE;
else if ( n == 2 ) // swap 3rd and 4th elements
{
{
int itmp;
itmp = src_add[2];
src_add[2] = src_add[3];
src_add[3] = itmp;
}
{
double dtmp[4];
for ( k = 0; k < 4; ++k ) dtmp[k] = weights[k][2];
for ( k = 0; k < 4; ++k ) weights[k][2] = weights[k][3];
for ( k = 0; k < 4; ++k ) weights[k][3] = dtmp[k];
}
}
}
for ( n = 0; n < 4; ++n )
{
rv->src_grid_add[nlink+n] = src_add[n];
rv->tgt_grid_add[nlink+n] = dst_add;
for ( k = 0; k < 4; ++k )
rv->wts[4*(nlink+n)+k] = weights[k][n];
}
} /* store_link_bicub */
/*
-----------------------------------------------------------------------
This routine computes the weights for a bicubic interpolation.
-----------------------------------------------------------------------
*/
void remap_bicub(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
{
/* Local variables */
int lwarn = TRUE;
int search_result;
long tgt_grid_size;
long n, icount;
long dst_add; /* destination addresss */
long iter; /* iteration counters */