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

remaplib: added get_srch_cells_reg2d()

parent 043eb230
Hide whitespace changes
Inline Side-by-side
src/interpol.c
View file @ a2225ff9
......@@ -372,32 +372,35 @@ void intconarr2(double missval, int lon_is_circular,
long i, ii = -1, jj = -1;
long gridsize1;
long nlon1 = nxm;
long nx, ny;
double findex = 0;
int *grid1_mask = NULL;
printf(" nxm, nym %ld %ld\n", nxm, nym);
nx = nxm - 1;
ny = nym - 1;
printf(" nx, ny %ld %ld\n", nx, ny);
//if ( lon_is_circular ) nlon1--;
gridsize1 = (nxm-1)*(nym-1);
gridsize1 = nx*ny;
deps = (int *) malloc(gridsize1*sizeof(int));
grid1_mask = (int *) calloc(1, gridsize1*sizeof(int));
for ( jj = 0; jj < nym-1; ++jj )
for ( ii = 0; ii < nxm-1; ++ii )
for ( jj = 0; jj < ny; ++jj )
for ( ii = 0; ii < nx; ++ii )
{
if ( !DBL_IS_EQUAL(fieldm[jj*(nxm-1)+ii], missval) ) grid1_mask[jj*(nxm-1)+ii] = 1;
if ( !DBL_IS_EQUAL(fieldm[jj*nx+ii], missval) ) grid1_mask[jj*nx+ii] = 1;
}
double grid1_bound_box[4];
grid1_bound_box[0] = ym[0];
grid1_bound_box[1] = ym[nym-1];
if ( ym[0] > ym[nym-1] )
grid1_bound_box[1] = ym[ny];
if ( ym[0] > ym[ny] )
{
grid1_bound_box[0] = ym[nym-1];
grid1_bound_box[0] = ym[ny];
grid1_bound_box[1] = ym[0];
}
grid1_bound_box[2] = xm[0];
grid1_bound_box[3] = xm[nxm-1];
grid1_bound_box[3] = xm[nx];
progressInit();
......@@ -479,6 +482,7 @@ void intconarr2(double missval, int lon_is_circular,
// for ( int k = 0; k < nxm; ++k ) printf("x: %d %g\n", k+1, xm[k]);
// for ( int k = 0; k < nym; ++k ) printf("y: %d %g\n", k+1, ym[k]*RAD2DEG);
long imin = nxm, imax = -1, jmin = nym, jmax = -1;
long im, jm;
lfound = rect_grid_search2(&jmin, &jmax, bound_box[0], bound_box[1], nym, ym);
bound_lon1 = bound_box[2];
......@@ -490,9 +494,9 @@ void intconarr2(double missval, int lon_is_circular,
if ( bound_lon2 > grid1_bound_box[3] && bound_lon1 < grid1_bound_box[3] ) bound_lon2 = grid1_bound_box[3];
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( long jm = jmin; jm <= jmax; ++jm )
for ( long im = imin; im <= imax; ++im )
deps[ndeps++] = jm*(nxm-1) + im;
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
deps[ndeps++] = jm*nx + im;
}
......@@ -507,9 +511,9 @@ void intconarr2(double missval, int lon_is_circular,
if ( bound_lon2 > grid1_bound_box[3] && bound_lon1 < grid1_bound_box[3] ) bound_lon2 = grid1_bound_box[3];
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( long jm = jmin; jm <= jmax; ++jm )
for ( long im = imin; im <= imax; ++im )
deps[ndeps++] = jm*(nxm-1) + im;
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
deps[ndeps++] = jm*nx + im;
}
bound_lon1 = bound_box[2];
......@@ -523,9 +527,9 @@ void intconarr2(double missval, int lon_is_circular,
if ( bound_lon2 > grid1_bound_box[3] && bound_lon1 < grid1_bound_box[3] ) bound_lon2 = grid1_bound_box[3];
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( long jm = jmin; jm <= jmax; ++jm )
for ( long im = imin; im <= imax; ++im )
deps[ndeps++] = jm*(nxm-1) + im;
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
deps[ndeps++] = jm*nx + im;
}
//for ( long id = 0; id < ndeps; ++id )
......@@ -566,8 +570,8 @@ void intconarr2(double missval, int lon_is_circular,
for ( int k = 0; k < num_deps; ++k )
{
int index1 = deps[k];
int ilat1 = index1/(nxm-1);
int ilon1 = index1 - ilat1*(nxm-1);
int ilat1 = index1/nx;
int ilon1 = index1 - ilat1*nx;
/*
if ( cdoVerbose )
printf(" dep: %d %d %d %d %d %d\n", k, nlonOut, nlatOut, index1, ilon1, ilat1);
......@@ -618,8 +622,8 @@ void intconarr2(double missval, int lon_is_circular,
{
int index1 = deps[k];
/*
int ilat1 = index1/(nxm-1);
int ilon1 = index1 - ilat1*(nxm-1);
int ilat1 = index1/nx;
int ilon1 = index1 - ilat1*nx;
long add1, add2;
add1 = index1;
......
src/remap.h
View file @ a2225ff9
......@@ -65,6 +65,9 @@ typedef struct {
double *reg2d_center_lon; /* reg2d lon/lat coordinates for */
double *reg2d_center_lat; /* each grid center in radians */
double *reg2d_corner_lon; /* reg2d lon/lat coordinates for */
double *reg2d_corner_lat; /* each grid corner in radians */
double *cell_center_lon; /* lon/lat coordinates for */
double *cell_center_lat; /* each grid center in radians */
double *cell_area; /* tot area of each grid cell */
......
src/remaplib.c
View file @ a2225ff9
......@@ -136,6 +136,8 @@ void remapGridFree(remapgrid_t *grid)
free(grid->mask);
if ( grid->reg2d_center_lat ) free(grid->reg2d_center_lat);
if ( grid->reg2d_center_lon ) free(grid->reg2d_center_lon);
if ( grid->reg2d_corner_lat ) free(grid->reg2d_corner_lat);
if ( grid->reg2d_corner_lon ) free(grid->reg2d_corner_lon);
free(grid->cell_center_lat);
free(grid->cell_center_lon);
if ( grid->cell_area ) free(grid->cell_area);
......@@ -330,6 +332,8 @@ void remapGridInitPointer(remapgrid_t *grid)
grid->mask = NULL;
grid->reg2d_center_lon = NULL;
grid->reg2d_center_lat = NULL;
grid->reg2d_corner_lon = NULL;
grid->reg2d_corner_lat = NULL;
grid->cell_center_lon = NULL;
grid->cell_center_lat = NULL;
grid->cell_area = NULL;
......@@ -792,17 +796,21 @@ void remap_reg2d_init(int gridID, remapgrid_t *grid)
{
char units[CDI_MAX_NAME];
long nx, nxm, ny;
long nxp1, nyp1;
nx = grid->dims[0];
ny = grid->dims[1];
nxp1 = nx+1;
nyp1 = ny+1;
nxm = nx;
if ( grid->is_cyclic ) nxm++;
if ( grid->size != nx*ny ) cdoAbort("Internal error, wrong dimensions!");
grid->reg2d_center_lon = (double *) realloc(grid->reg2d_center_lon, (nx+1)*sizeof(double));
grid->reg2d_center_lat = (double *) realloc(grid->reg2d_center_lat, ny*sizeof(double));
grid->reg2d_center_lon = (double *) realloc(grid->reg2d_center_lon, nxm*sizeof(double));
grid->reg2d_center_lat = (double *) realloc(grid->reg2d_center_lat, ny*sizeof(double));
gridInqXvals(gridID, grid->reg2d_center_lon);
gridInqYvals(gridID, grid->reg2d_center_lat);
......@@ -991,7 +999,7 @@ void remap_grids_init(int map_type, int lextrapolate, int gridID1, remapgrid_t *
src_grid->remap_grid_type = -1;
tgt_grid->remap_grid_type = -1;
if ( (map_type == MAP_TYPE_BILINEAR || map_type == MAP_TYPE_BICUBIC || map_type == MAP_TYPE_DISTWGT ) &&
if ( (map_type == MAP_TYPE_BILINEAR || map_type == MAP_TYPE_BICUBIC || map_type == MAP_TYPE_DISTWGT || map_type == MAP_TYPE_CONTEST ) &&
!gridIsRotated(gridID1) &&
(gridInqType(gridID1) == GRID_LONLAT || gridInqType(gridID1) == GRID_GAUSSIAN) )
src_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
......@@ -4264,9 +4272,107 @@ void store_link_cnsrv(remapvars_t *rv, long add1, long add2, double *restrict we
} /* store_link_cnsrv */
static
long get_srch_cells_reg2d(const int *restrict src_grid_dims,
const double *restrict src_corner_lat, const double *restrict src_corner_lon,
const restr_t *tgt_cell_bound_box, int *srch_add)
{
long num_srch_cells; /* num cells in restricted search arrays */
long min_add; /* addresses for restricting search of */
long max_add; /* destination grid */
long n, n2; /* generic counters */
long src_grid_add; /* current linear address for src cell */
long tgt_grid_addm4, src_grid_addm4;
int lmask;
long nx, ny;
long nxp1, nyp1;
double grid1_bound_box[4];
double src_lon_min, src_lon_max;
nx = src_grid_dims[0];
ny = src_grid_dims[1];
nxp1 = nx+1;
nyp1 = ny+1;
/*
grid1_bound_box[0] = src_corner_lat[0];
grid1_bound_box[1] = src_corner_lat[ny];
if ( ym[0] > ym[ny] )
{
grid1_bound_box[0] = src_corner_lat[ny];
grid1_bound_box[1] = src_corner_lat[0];
}
grid1_bound_box[2] = src_corner_lon[0];
grid1_bound_box[3] = src_corner_lon[nx];
src_lon_min = src_corner_lon[0];
src_lon_max = src_corner_lon[nx];
double bound_lon1, bound_lon2;
double bound_box[4];
boundbox_from_corners(i, nc2, x, y, bound_box);
restrict_boundbox(grid1_bound_box, bound_box);
bound_lon1 = bound_box[2];
bound_lon2 = bound_box[3];
*/
num_srch_cells = 0;
/*
long imin = nxm, imax = -1, jmin = nym, jmax = -1;
long im, jm;
lfound = rect_grid_search2(&jmin, &jmax, bound_box[0], bound_box[1], nym, ym);
bound_lon1 = bound_box[2];
bound_lon2 = bound_box[3];
if ( bound_lon1 <= src_lon_max && bound_lon2 >= src_lon_min )
{
//printf("b1 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
srch_add[num_srch_cells++] = jm*nx + im;
}
bound_lon1 = bound_box[2];
bound_lon2 = bound_box[3];
if ( bound_lon1 <= src_lon_min && bound_lon2 > src_lon_min )
{
bound_lon1 += 2*M_PI;
bound_lon2 += 2*M_PI;
//printf("b2 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
srch_add[num_srch_cells++] = jm*nx + im;
}
bound_lon1 = bound_box[2];
bound_lon2 = bound_box[3];
if ( bound_lon1 < src_lon_max && bound_lon2 >= src_lon_max )
{
bound_lon1 -= 2*M_PI;
bound_lon2 -= 2*M_PI;
//printf("b3 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
srch_add[num_srch_cells++] = jm*nx + im;
}
*/
return (num_srch_cells);
}
static
long get_srch_cells(long tgt_grid_add, long nbins, int *bin_addr1, int *bin_addr2,
restr_t *tgt_cell_bound_box, restr_t *src_cell_bound_box, long grid2_size, int *srch_add)
restr_t *tgt_cell_bound_box, restr_t *src_cell_bound_box, long src_grid_size, int *srch_add)
{
long num_srch_cells; /* num cells in restricted search arrays */
long min_add; /* addresses for restricting search of */
......@@ -4279,7 +4385,7 @@ long get_srch_cells(long tgt_grid_add, long nbins, int *bin_addr1, int *bin_addr
/* Restrict searches first using search bins */
min_add = grid2_size - 1;
min_add = src_grid_size - 1;
max_add = 0;
for ( n = 0; n < nbins; ++n )
......@@ -4525,7 +4631,7 @@ void remap_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
/* Get search cells */
num_srch_cells = get_srch_cells(src_grid_add, nbins, src_grid->bin_addr, tgt_grid->bin_addr,
src_grid->cell_bound_box, tgt_grid->cell_bound_box , grid2_size, srch_add);
src_grid->cell_bound_box, tgt_grid->cell_bound_box, grid2_size, srch_add);
if ( num_srch_cells == 0 ) continue;
......@@ -4751,7 +4857,7 @@ void remap_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
/* Get search cells */
num_srch_cells = get_srch_cells(tgt_grid_add, nbins, tgt_grid->bin_addr, src_grid->bin_addr,
tgt_grid->cell_bound_box, src_grid->cell_bound_box , grid1_size, srch_add);
tgt_grid->cell_bound_box, src_grid->cell_bound_box, grid1_size, srch_add);
if ( num_srch_cells == 0 ) continue;
......@@ -5328,6 +5434,7 @@ void remap_contest(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
grid_store_t *grid_store = NULL;
double findex = 0;
long num_weights = 0;
int remap_grid_type = src_grid->remap_grid_type;
progressInit();
......@@ -5455,8 +5562,16 @@ void remap_contest(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
/* Get search cells */
num_srch_cells = get_srch_cells(tgt_grid_add, nbins, tgt_grid->bin_addr, src_grid->bin_addr,
tgt_grid->cell_bound_box, src_grid->cell_bound_box , grid1_size, srch_add);
/*
if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
num_srch_cells = get_srch_cells_reg2d(src_grid->dims, src_grid->reg2d_corner_lat, src_grid->reg2d_corner_lon,
tgt_cell_bound_box, srch_add);
else
*/
num_srch_cells = get_srch_cells(tgt_grid_add, nbins, tgt_grid->bin_addr, src_grid->bin_addr,
tgt_grid->cell_bound_box, src_grid->cell_bound_box, grid1_size, srch_add);
printf("tgt_grid_add %ld num_srch_cells %ld\n", tgt_grid_add, num_srch_cells);
if ( num_srch_cells == 0 ) continue;
......@@ -5549,7 +5664,7 @@ void remap_contest(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
}
}
//correct_weights(num_weights, weight);
correct_weights(num_weights, weight);
#endif
for ( n = 0; n < num_weights; ++n )
......
Supports Markdown
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

Imprint | Privacy Policy