Commit 9f31bcd9 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remapdis/remapnn without extrapolation on non global curvilinear grids:...

remapdis/remapnn without extrapolation on non global curvilinear grids: replaced expansion of borders by local search.
parent f9c9977a
2017-12-06 Uwe Schulzweida
* remapdis/remapnn without extrapolation on non global curvilinear grids:
replaced expansion of borders by local search
2017-11-30 Uwe Schulzweida
* New operator not: logical NOT (1, if x equal 0; else 0)
......
......@@ -434,19 +434,6 @@ void get_remap_env(void)
}
}
static
void set_halo_to_missval(size_t nx, size_t ny, double *array, double missval)
{
for ( size_t j = 0; j < ny+4; j++ ) array[j*(nx+4)+0] = missval;
for ( size_t j = 0; j < ny+4; j++ ) array[j*(nx+4)+1] = missval;
for ( size_t j = 0; j < ny+4; j++ ) array[j*(nx+4)+nx+2] = missval;
for ( size_t j = 0; j < ny+4; j++ ) array[j*(nx+4)+nx+3] = missval;
for ( size_t i = 0; i < nx+4; i++ ) array[ 0*(nx+4)+i] = missval;
for ( size_t i = 0; i < nx+4; i++ ) array[ 1*(nx+4)+i] = missval;
for ( size_t i = 0; i < nx+4; i++ ) array[(ny+2)*(nx+4)+i] = missval;
for ( size_t i = 0; i < nx+4; i++ ) array[(ny+3)*(nx+4)+i] = missval;
}
static
bool is_global_grid(int gridID)
{
......@@ -901,12 +888,6 @@ void *Remap(void *argument)
if ( map_type == MAP_TYPE_DISTWGT && !lextrapolate ) remap_extrapolate = true;
if ( gridIsCircular(gridID1) && !lextrapolate ) remap_extrapolate = true;
if ( map_type == MAP_TYPE_DISTWGT && !remap_extrapolate && gridInqSize(gridID1) > 1 && !is_global_grid(gridID1) )
{
remaps[0].gridsize += 4*(gridInqXsize(gridID1)+2) + 4*(gridInqYsize(gridID1)+2);
remaps[0].src_grid.non_global = true;
}
if ( gridInqType(gridID1) == GRID_GME ) gridsize = remaps[0].src_grid.nvgp;
if ( gridsize != remaps[0].gridsize )
......@@ -1016,28 +997,6 @@ void *Remap(void *argument)
gridsize = gridInqSize(gridID1);
if ( gridIsCircular(gridID1) && !lextrapolate ) remap_extrapolate = true;
if ( map_type == MAP_TYPE_DISTWGT && !remap_extrapolate && gridInqSize(gridID1) > 1 && !is_global_grid(gridID1) )
{
if ( cdoVerbose ) cdoPrint("---> Expand array!");
long nx = gridInqXsize(gridID1);
long ny = gridInqYsize(gridID1);
size_t gridsize_new = gridsize + 4*(nx+2) + 4*(ny+2);
if ( gridsize_new > grid1sizemax )
{
grid1sizemax = gridsize_new;
array1 = (double*) Realloc(array1, grid1sizemax*sizeof(double));
imask = (int*) Realloc(imask, grid1sizemax*sizeof(int));
}
for ( long j = ny-1; j >= 0; j-- )
for ( long i = nx-1; i >= 0; i-- )
array1[(j+2)*(nx+4)+i+2] = array1[j*nx+i];
set_halo_to_missval(nx, ny, array1, missval);
gridsize = gridsize_new;
nmiss1 += 4*(nx+2) + 4*(ny+2);
}
for ( size_t i = 0; i < gridsize; i++ )
imask[i] = DBL_IS_EQUAL(array1[i], missval) ? FALSE : TRUE;
......
......@@ -322,173 +322,6 @@ void *gs_create_full(size_t n, const double *restrict lons, const double *restri
return (void*) full;
}
static
void cal_bound_box(struct gridsearch *gs, size_t n, const double *restrict lons, const double *restrict lats)
{
enum {mlon=180, mlat=90};
const size_t masksize = mlon*mlat;
bool *mask1 = (bool*) Malloc(masksize*sizeof(bool));
for ( size_t i = 0; i < masksize; ++i ) mask1[i] = 0;
for ( size_t i = 0; i < n; ++i )
{
double fi = lons[i];
double fj = lats[i]+0.5*M_PI;
if ( fi < 0 ) fi += 2*M_PI;
fi *= mlon/(2*M_PI);
fj *= mlat/M_PI;
int64_t ii = (int64_t) llround(fi);
int64_t jj = (int64_t) llround(fj);
if ( i < 1000 )
printf("lon=%g/lat=%g ii=%ld jj=%ld fi=%g fj=%g\n", lons[i]*RAD2DEG, lats[i]*RAD2DEG, (long)ii, (long)jj, fi, fj);
mask1[jj*mlon+ii] = 1;
}
double xlonmin = 999., xlonmax = 999.;
double xlatmin = 999., xlatmax = 999.;
for ( size_t j = 0; j < mlat; ++j )
{
for ( size_t i = 0; i < mlon; ++i )
if ( mask1[j*mlon+i] )
{
xlatmin = (j+0.)*M_PI/mlat - 0.5*M_PI;
break;
}
if ( IS_NOT_EQUAL(xlatmin, 999.) ) break;
}
for ( size_t jj = 0; jj < mlat; ++jj )
{
size_t j = mlat-jj-1;
for ( size_t i = 0; i < mlon; ++i )
if ( mask1[j*mlon+i] )
{
xlatmax = (j+1.)*M_PI/mlat - 0.5*M_PI;
break;
}
if ( IS_NOT_EQUAL(xlatmax, 999.) ) break;
}
size_t ii = mlon;
for ( size_t i = 0; i < mlon; ++i )
{
for ( size_t j = 0; j < mlat; ++j )
if ( mask1[j*mlon+i] )
{
ii = i;
xlonmin = (i+.0)*2*M_PI/mlon;
break;
}
if ( IS_NOT_EQUAL(xlonmin, 999.) ) break;
}
if ( cdoVerbose ) printf("ii= %zu\n", ii);
for ( size_t i = ii; i < mlon; ++i )
{
ii = mlon;
size_t imask = 0;
for ( size_t j = 0; j < mlat; ++j )
{
if ( mask1[j*mlon+i] )
{
imask++;
xlonmax = (i+1.)*2*M_PI/mlon;
}
}
if ( imask == 0 )
{
ii = i+1;
break;
}
}
if ( cdoVerbose ) cdoPrint("boundbox: lonmin=%g lonmax=%g latmin=%g latmax=%g",
xlonmin*RAD2DEG, xlonmax*RAD2DEG, xlatmin*RAD2DEG, xlatmax*RAD2DEG);
if ( cdoVerbose ) printf("ii= %zu\n", ii);
for ( size_t i = ii; i < mlon; ++i )
{
bool lfound = false;
for ( size_t j = 0; j < mlat; ++j )
if ( mask1[j*mlon+i] )
{
lfound = true;
xlonmin = (i+.0)*2*M_PI/mlon;
break;
}
if ( lfound ) break;
}
if ( xlonmin > xlonmax ) xlonmin -= 2*M_PI;
if ( cdoVerbose ) cdoPrint("boundbox: lonmin=%g lonmax=%g latmin=%g latmax=%g",
xlonmin*RAD2DEG, xlonmax*RAD2DEG, xlatmin*RAD2DEG, xlatmax*RAD2DEG);
for ( size_t j = 0; j < mlat; ++j )
{
for ( size_t i = 0; i < mlon; ++i )
{
printf("%1d", mask1[j*mlon+i]);
}
printf("\n");
}
double lonmin = 1.e33, lonmax = -1.e33;
double latmin = 1.e33, latmax = -1.e33;
for ( size_t i = 0; i < n; ++i )
{
if ( lons[i] < lonmin ) lonmin = lons[i];
if ( lons[i] > lonmax ) lonmax = lons[i];
if ( lats[i] < latmin ) latmin = lats[i];
if ( lats[i] > latmax ) latmax = lats[i];
}
gs->lonmin = xlonmin;
gs->lonmax = xlonmax;
gs->latmin = xlatmin;
gs->latmax = xlatmax;
if ( cdoVerbose ) cdoPrint("boundbox: lonmin=%g lonmax=%g latmin=%g latmax=%g",
lonmin*RAD2DEG, lonmax*RAD2DEG, latmin*RAD2DEG, latmax*RAD2DEG);
if ( mask1 ) free(mask1);
}
static
void cal_mask(struct gridsearch *gs)
{
enum {mmin = 100};
const double fact = 1;
double dlon = gs->lonmax - gs->lonmin;
double dlat = gs->latmax - gs->latmin;
if ( cdoVerbose ) cdoPrint("mask: dlon=%g, dlat=%g", dlon*RAD2DEG, dlat*RAD2DEG);
size_t sqrtn = (size_t) sqrt((double)gs->n);
if ( cdoVerbose ) cdoPrint("n=%zu sqrt(n)=%zu", gs->n, sqrtn);
size_t mlat = fact*(sqrtn/2.)*(1+dlat/dlon);
size_t mlon = mlat*dlon/dlat;
if ( cdoVerbose ) cdoPrint("mlon=%zu mlat=%zu mn=%zu", mlon, mlat, mlon*mlat);
}
void gridsearch_bound_poly(struct gridsearch *gs, size_t dims[2], size_t n, const double *restrict lons, const double *restrict lats)
{
printf("%s\n", __func__);
size_t nx = dims[0];
size_t ny = dims[1];
printf("nx %zu ny %zu\n", nx, ny);
if ( n == 1 || nx*ny != n ) return;
#define PMAX 40
size_t pmax = PMAX;
size_t pn = 0;
double px[PMAX];
double py[PMAX];
for ( size_t i = 0; i < nx; ++i )
printf("x1 %zu lon/lat %g/%g\n", i+1, lons[0*nx+i]*RAD2DEG, lats[0*nx+i]*RAD2DEG);
for ( size_t i = 0; i < nx; ++i )
printf("x2 %zu lon/lat %g/%g\n", i+1, lons[(ny-1)*nx+i]*RAD2DEG, lats[(ny-1)*nx+i]*RAD2DEG);
for ( size_t i = 0; i < ny; ++i )
printf("y1 %zu lon/lat %g/%g\n", i+1, lons[i*nx+0]*RAD2DEG, lats[i*nx+0]*RAD2DEG);
for ( size_t i = 0; i < ny; ++i )
printf("y2 %zu lon/lat %g/%g\n", i+1, lons[i*nx+(nx-1)]*RAD2DEG, lats[i*nx+(nx-1)]*RAD2DEG);
}
struct gridsearch *gridsearch_create(bool xIsCyclic, size_t dims[2], size_t n, const double *restrict lons, const double *restrict lats)
{
......@@ -503,6 +336,9 @@ struct gridsearch *gridsearch_create(bool xIsCyclic, size_t dims[2], size_t n, c
gs->n = n;
if ( n == 0 ) return gs;
gs->plons = lons;
gs->plats = lats;
if ( gs->method_nn == GridsearchMethod::kdtree ) gs->search_container = gs_create_kdtree(n, lons, lats, gs);
else if ( gs->method_nn == GridsearchMethod::nanoflann ) gs->search_container = gs_create_nanoflann(n, lons, lats, gs);
......@@ -534,9 +370,6 @@ struct gridsearch *gridsearch_create_nn(bool xIsCyclic, size_t dims[2], size_t n
gs->search_radius = cdo_default_search_radius();
// cal_bound_box(gs, n, lons, lats);
// cal_mask(gs);
return gs;
}
......@@ -606,12 +439,42 @@ size_t gs_nearest_kdtree(void *search_container, double lon, double lat, double
return index;
}
//#define TEST_QUAD
#ifdef TEST_QUAD
bool point_in_quad(bool is_cyclic, size_t nx, size_t ny, size_t i, size_t j, size_t adds[4], double lons[4], double lats[4],
double plon, double plat, const double *restrict center_lon, const double *restrict center_lat);
#endif
static
size_t llindex_in_quad(struct gridsearch *gs, size_t index, double lon, double lat)
{
size_t ret_index = index;
if ( ret_index != GS_NOT_FOUND )
{
ret_index = GS_NOT_FOUND;
size_t nx = gs->dims[0];
size_t ny = gs->dims[1];
size_t adds[4];
double lons[4];
double lats[4];
bool is_cyclic = gs->is_cyclic;
for ( unsigned k = 0; k < 4; ++k )
{
/* Determine neighbor addresses */
size_t j = index/nx;
size_t i = index - 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, adds, lons, lats,
lon, lat, gs->plons, gs->plats) )
{
ret_index = index;
break;
}
}
}
return ret_index;
}
static
size_t gs_nearest_nanoflann(void *search_container, double lon, double lat, double *prange, struct gridsearch *gs)
......@@ -638,41 +501,8 @@ size_t gs_nearest_nanoflann(void *search_container, double lon, double lat, doub
nft->findNeighbors(resultSet, query_pt, nanoflann::SearchParams(10));
//printf("%zu %g\n", ret_index, out_dist_sqr);
#ifdef TEST_QUAD
if ( !gs->extrapolate && gs->is_curve )
{
if ( ret_index != GS_NOT_FOUND )
{
size_t nx = gs->dims[0];
size_t ny = gs->dims[1];
size_t adds[4];
double lons[4];
double lats[4];
bool is_cyclic = gs->is_cyclic;
for ( unsigned k = 0; k < 4; ++k )
{
/* Determine neighbor addresses */
size_t j = ret_index/nx;
size_t i = ret_index - 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, adds, lons, lats,
lon, lat, gs->plons, gs->plats) )
{
index = ret_index;
break;
}
}
}
}
else
#else
{
index = ret_index;
*prange = out_dist_sqr;
}
#endif
index = ret_index;
*prange = out_dist_sqr;
//float frange = range;
//if ( !(frange < range0) ) node = NULL;
......@@ -735,6 +565,8 @@ size_t gridsearch_nearest(struct gridsearch *gs, double lon, double lat, double
else if ( gs->method_nn == GridsearchMethod::full ) index = gs_nearest_full(sc, lon, lat, prange);
else cdoAbort("%s::method_nn undefined!", __func__);
// clang-format on
if ( !gs->extrapolate && gs->is_curve ) index = llindex_in_quad(gs, index, lon, lat);
}
return index;
......@@ -748,13 +580,17 @@ size_t gs_qnearest_kdtree(struct gridsearch *gs, double lon, double lat, double
kdTree_t *kdt = (kdTree_t *) gs->search_container;
if ( kdt == NULL ) return nadds;
kdata_t query_pt[3];
kdata_t range0 = gs_set_range(prange);
kdata_t range = range0;
struct pqueue *result = NULL;
kdata_t query_pt[3];
LLtoXYZ<kdata_t>(lon, lat, query_pt);
if ( !gs->extrapolate )
for ( unsigned j = 0; j < 3; ++j )
if ( query_pt[j] < gs->min[j] || query_pt[j] > gs->max[j] ) return nadds;
if ( gs )
{
result = kd_qnearest(kdt->node, query_pt, &range, nnn, 3);
......@@ -803,6 +639,10 @@ size_t gs_qnearest_nanoflann(struct gridsearch *gs, double lon, double lat, doub
NFDATATYPE query_pt[3];
LLtoXYZ<NFDATATYPE>(lon, lat, query_pt);
if ( !gs->extrapolate )
for ( unsigned j = 0; j < 3; ++j )
if ( query_pt[j] < gs->min[j] || query_pt[j] > gs->max[j] ) return nadds;
if ( gs )
{
std::vector<NFDATATYPE> out_dist_sqr(nnn);
......@@ -829,6 +669,23 @@ size_t gridsearch_qnearest(struct gridsearch *gs, double lon, double lat, double
else if ( gs->method_nn == GridsearchMethod::nanoflann ) nadds = gs_qnearest_nanoflann(gs, lon, lat, prange, nnn, adds, dist);
else cdoAbort("%s::method_nn undefined!", __func__);
// clang-format on
if ( !gs->extrapolate && gs->is_curve )
{
size_t nadds_old = nadds;
nadds = 0;
for ( size_t i = 0; i < nadds_old; ++i )
{
size_t index = adds[i];
index = llindex_in_quad(gs, index, lon, lat);
if ( index != GS_NOT_FOUND )
{
dist[nadds] = dist[i];
adds[nadds] = adds[i];
nadds++;
}
}
}
}
return nadds;
......
......@@ -450,71 +450,6 @@ int expand_lonlat_grid(int gridID)
return gridIDnew;
}
static
int expand_curvilinear_grid(int gridID)
{
if ( cdoVerbose ) cdoPrint("expand_curvilinear_grid");
size_t gridsize = gridInqSize(gridID);
long nx = (long) gridInqXsize(gridID);
long ny = (long) gridInqYsize(gridID);
long nxp4 = nx+4;
long nyp4 = ny+4;
size_t gridsize_new = gridsize + 4*(nx+2) + 4*(ny+2);
double *xvals = (double*) Malloc(gridsize_new*sizeof(double));
double *yvals = (double*) Malloc(gridsize_new*sizeof(double));
gridInqXvals(gridID, xvals);
gridInqYvals(gridID, yvals);
int gridIDnew = gridCreate(GRID_CURVILINEAR, nxp4*nyp4);
gridDefXsize(gridIDnew, nxp4);
gridDefYsize(gridIDnew, nyp4);
grid_copy_attributes(gridID, gridIDnew);
for ( long j = ny-1; j >= 0; j-- )
for ( long i = nx-1; i >= 0; i-- )
xvals[(j+2)*(nx+4)+i+2] = xvals[j*nx+i];
for ( long j = ny-1; j >= 0; j-- )
for ( long i = nx-1; i >= 0; i-- )
yvals[(j+2)*(nx+4)+i+2] = yvals[j*nx+i];
for ( long j = 2; j < nyp4-2; j++ )
{
xvals[j*nxp4 ] = intlin(3.0, xvals[j*nxp4+3], 0.0, xvals[j*nxp4+2], 1.0);
xvals[j*nxp4+1] = intlin(2.0, xvals[j*nxp4+3], 0.0, xvals[j*nxp4+2], 1.0);
yvals[j*nxp4 ] = intlin(3.0, yvals[j*nxp4+3], 0.0, yvals[j*nxp4+2], 1.0);
yvals[j*nxp4+1] = intlin(2.0, yvals[j*nxp4+3], 0.0, yvals[j*nxp4+2], 1.0);
xvals[j*nxp4+nxp4-2] = intlin(2.0, xvals[j*nxp4+nxp4-4], 0.0, xvals[j*nxp4+nxp4-3], 1.0);
xvals[j*nxp4+nxp4-1] = intlin(3.0, xvals[j*nxp4+nxp4-4], 0.0, xvals[j*nxp4+nxp4-3], 1.0);
yvals[j*nxp4+nxp4-2] = intlin(2.0, yvals[j*nxp4+nxp4-4], 0.0, yvals[j*nxp4+nxp4-3], 1.0);
yvals[j*nxp4+nxp4-1] = intlin(3.0, yvals[j*nxp4+nxp4-4], 0.0, yvals[j*nxp4+nxp4-3], 1.0);
}
for ( long i = 0; i < nxp4; i++ )
{
xvals[0*nxp4+i] = intlin(3.0, xvals[3*nxp4+i], 0.0, xvals[2*nxp4+i], 1.0);
xvals[1*nxp4+i] = intlin(2.0, xvals[3*nxp4+i], 0.0, xvals[2*nxp4+i], 1.0);
yvals[0*nxp4+i] = intlin(3.0, yvals[3*nxp4+i], 0.0, yvals[2*nxp4+i], 1.0);
yvals[1*nxp4+i] = intlin(2.0, yvals[3*nxp4+i], 0.0, yvals[2*nxp4+i], 1.0);
xvals[(nyp4-2)*nxp4+i] = intlin(2.0, xvals[(nyp4-4)*nxp4+i], 0.0, xvals[(nyp4-3)*nxp4+i], 1.0);
xvals[(nyp4-1)*nxp4+i] = intlin(3.0, xvals[(nyp4-4)*nxp4+i], 0.0, xvals[(nyp4-3)*nxp4+i], 1.0);
yvals[(nyp4-2)*nxp4+i] = intlin(2.0, yvals[(nyp4-4)*nxp4+i], 0.0, yvals[(nyp4-3)*nxp4+i], 1.0);
yvals[(nyp4-1)*nxp4+i] = intlin(3.0, yvals[(nyp4-4)*nxp4+i], 0.0, yvals[(nyp4-3)*nxp4+i], 1.0);
}
gridDefXvals(gridIDnew, xvals);
gridDefYvals(gridIDnew, yvals);
Free(xvals);
Free(yvals);
return gridIDnew;
}
/*****************************************************************************/
static
......@@ -847,13 +782,6 @@ void remap_grids_init(int map_type, bool lextrapolate, int gridID1, remapgrid_t
src_grid->gridID = gridID1 = gridToCurvilinear(src_grid->gridID, lbounds);
}
if ( !src_grid->lextrapolate && gridInqSize(src_grid->gridID) > 1 &&
map_type == MAP_TYPE_DISTWGT &&
(gridInqType(gridID1) == GRID_CURVILINEAR && src_grid->non_global) )
{
src_grid->gridID = gridID1 = expand_curvilinear_grid(gridID1);
}
//if ( src_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
remap_define_grid(map_type, gridID1, src_grid, "Source");
remap_define_grid(map_type, gridID2, tgt_grid, "Target");
......
......@@ -5,23 +5,6 @@ CDO=cdo
FORMAT="-f srv -b F32"
########################################################################
#
# Remap
#
cdo -f grb setrtomiss,0,10000 -gridboxmean,8,8 -topo bathy4.grb
#
GRIDS="n16 n32"
RMODS="bil bic dis nn con con2 laf"
RMODS="con2"
IFILE=bathy4.grb
for GRID in $GRIDS; do
for RMOD in $RMODS; do
OFILE=${GRID}_${RMOD}
$CDO $FORMAT remap${RMOD},$GRID $IFILE ${OFILE}_ref
done
done
exit
########################################################################
#
# Remap regional grid
#
GRID=spain.grid
......@@ -47,6 +30,23 @@ done
exit
########################################################################
#
# Remap
#
cdo -f grb setrtomiss,0,10000 -gridboxmean,8,8 -topo bathy4.grb
#
GRIDS="n16 n32"
RMODS="bil bic dis nn con con2 laf"
RMODS="con2"
IFILE=bathy4.grb
for GRID in $GRIDS; do
for RMOD in $RMODS; do
OFILE=${GRID}_${RMOD}
$CDO $FORMAT remap${RMOD},$GRID $IFILE ${OFILE}_ref
done
done
exit
########################################################################
#
# Test GRIB files
#
IFILES="testfile01 testfile02 testfile03"
......
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