Commit e09eb3fd authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

New operator: setmisstodis - Set missing value to the distance-weighted...

New operator: setmisstodis - Set missing value to the distance-weighted average of the nearest neighbors
parent b024cebb
......@@ -3,6 +3,10 @@
* using CDI library version 1.7.0
* Version 1.7.0 released
2015-10-23 Uwe Schulzweida
* New operator: setmisstodis - Set missing value to the distance-weighted average of the nearest neighbors
2015-10-17 Uwe Schulzweida
* diff: set checkrel=TRUE
......
......@@ -14,6 +14,7 @@ Version 1.7.0 (28 October 2015):
* remapycon: First order conservative remapping (new implementation of remapcon)
* genycon: Generate 1st order conservative remap weights (new implementation of gencon)
* setmisstonn: Set missing value to nearest neightbor
* setmisstodis: Set missing value to the distance-weighted average of the nearest neighbors
* ap2pl: Interpolate 3D variables on hybrid sigma height coordinates to pressure levels
* vertstd1: Vertical standard deviation [Divisor is (n-1)]
* vertvar1: Vertical variance [Divisor is (n-1)]
......
......@@ -271,15 +271,18 @@ void fillmiss_one_step(field_t *field1, field_t *field2, int maxfill)
}
void setmisstonn(field_t *field1, field_t *field2, int maxfill)
void grid_search_nbr(struct gridsearch *gs, int num_neighbors, int *restrict nbr_add, double *restrict nbr_dist, double plon, double plat);
double nbr_compute_weights(unsigned num_neighbors, const int *restrict src_grid_mask, int *restrict nbr_mask, const int *restrict nbr_add, double *restrict nbr_dist);
unsigned nbr_normalize_weights(unsigned num_neighbors, double dist_tot, const int *restrict nbr_mask, int *restrict nbr_add, double *restrict nbr_dist);
static
void setmisstodis(field_t *field1, field_t *field2, int num_neighbors)
{
int gridID = field1->grid;
double missval = field1->missval;
double *array1 = field1->ptr;
double *array2 = field2->ptr;
UNUSED(maxfill);
unsigned gridsize = gridInqSize(gridID);
unsigned nmiss = field1->nmiss;
......@@ -319,6 +322,7 @@ void setmisstonn(field_t *field1, field_t *field2, int maxfill)
unsigned nv = 0, nm = 0;
for ( unsigned i = 0; i < gridsize; ++i )
{
array2[i] = array1[i];
if ( DBL_IS_EQUAL(array1[i], missval) )
{
mindex[nm] = i;
......@@ -326,7 +330,6 @@ void setmisstonn(field_t *field1, field_t *field2, int maxfill)
}
else
{
array2[i] = array1[i];
lons[nv] = xvals[i];
lats[nv] = yvals[i];
vindex[nv] = i;
......@@ -335,14 +338,23 @@ void setmisstonn(field_t *field1, field_t *field2, int maxfill)
}
if ( nv != nvals ) cdoAbort("Internal problem, number of valid values differ!");
int nbr_mask[num_neighbors]; /* mask at nearest neighbors */
int nbr_add[num_neighbors]; /* source address at nearest neighbors */
double nbr_dist[num_neighbors]; /* angular distance four nearest neighbors */
clock_t start, finish;
start = clock();
struct gridsearch *gs = NULL;
if ( nmiss ) gs = gridsearch_create_nn(nvals, lons, lats);
if ( nmiss )
{
if ( num_neighbors == 1 )
gs = gridsearch_create_nn(nvals, lons, lats);
else
gs = gridsearch_create(nvals, lons, lats);
}
finish = clock();
......@@ -352,12 +364,10 @@ void setmisstonn(field_t *field1, field_t *field2, int maxfill)
start = clock();
extern double gridsearch_radius;
double findex = 0;
double search_radius = gridsearch_radius*DEG2RAD;
double range = SQR(2*search_radius);
#pragma omp parallel for default(none) shared(findex, mindex, vindex, array1, array2, xvals, yvals, range, gs, nmiss)
#pragma omp parallel for default(none) shared(findex, mindex, vindex, array1, array2, xvals, yvals, gs, nmiss, num_neighbors) \
private(nbr_mask, nbr_add, nbr_dist)
for ( unsigned i = 0; i < nmiss; ++i )
{
#if defined(_OPENMP)
......@@ -366,14 +376,19 @@ void setmisstonn(field_t *field1, field_t *field2, int maxfill)
findex++;
if ( cdo_omp_get_thread_num() == 0 ) progressStatus(0, 1, findex/nmiss);
double prange = range;
unsigned index = gridsearch_nearest(gs, xvals[mindex[i]], yvals[mindex[i]], &prange);
if ( index == GS_NOT_FOUND ) index = mindex[i];
else index = vindex[index];
grid_search_nbr(gs, num_neighbors, nbr_add, nbr_dist, xvals[mindex[i]], yvals[mindex[i]]);
// printf("%u %u %d\n", i, index, (int)(prange*100000));
/* Compute weights based on inverse distance if mask is false, eliminate those points */
double dist_tot = nbr_compute_weights(num_neighbors, NULL, nbr_mask, nbr_add, nbr_dist);
array2[mindex[i]] = array1[index];
/* Normalize weights and store the link */
unsigned nadds = nbr_normalize_weights(num_neighbors, dist_tot, nbr_mask, nbr_add, nbr_dist);
if ( nadds )
{
double result = 0;
for ( unsigned n = 0; n < nadds; ++n ) result += array1[vindex[nbr_add[n]]]*nbr_dist[n];
array2[mindex[i]] = result;
}
}
if ( mindex ) Free(mindex);
......@@ -394,20 +409,19 @@ void setmisstonn(field_t *field1, field_t *field2, int maxfill)
void *Fillmiss(void *argument)
{
int gridID;
int nrecs;
int recID, varID, levelID;
int nmiss, i, nfill = 1;
int nrecs, recID, varID, levelID;
void (*fill_method) (field_t *fin , field_t *fout , int) = NULL;
cdoInitialize(argument);
int FILLMISS = cdoOperatorAdd("fillmiss" , 0, 0, "nfill");
int FILLMISSONESTEP = cdoOperatorAdd("fillmiss2" , 0, 0, "nfill");
int SETMISSTONN = cdoOperatorAdd("setmisstonn" , 0, 0, "nfill");
int FILLMISS = cdoOperatorAdd("fillmiss" , 0, 0, "nfill");
int FILLMISSONESTEP = cdoOperatorAdd("fillmiss2" , 0, 0, "nfill");
int SETMISSTONN = cdoOperatorAdd("setmisstonn" , 0, 0, "");
int SETMISSTODIS = cdoOperatorAdd("setmisstodis" , 0, 0, "number of neighbors");
int operatorID = cdoOperatorID();
int nfill = 1;
if ( operatorID == FILLMISS )
{
fill_method = &fillmiss;
......@@ -418,7 +432,12 @@ void *Fillmiss(void *argument)
}
else if ( operatorID == SETMISSTONN )
{
fill_method = &setmisstonn;
fill_method = &setmisstodis;
}
else if ( operatorID == SETMISSTODIS )
{
nfill = 4;
fill_method = &setmisstodis;
}
/* Argument handling */
......@@ -428,8 +447,8 @@ void *Fillmiss(void *argument)
if ( oargc == 1 )
{
nfill = atoi(oargv[0]);
if ( operatorID != FILLMISSONESTEP )
nfill = parameter2int(oargv[0]);
if ( operatorID == FILLMISS )
{
if ( nfill < 1 || nfill > 4 ) cdoAbort("nfill out of range!");
}
......@@ -479,11 +498,10 @@ void *Fillmiss(void *argument)
}
else
{
gridID = vlistInqVarGrid(vlistID1, varID);
int gridID = vlistInqVarGrid(vlistID1, varID);
if ( operatorID != SETMISSTONN &&
(gridInqType(gridID) == GRID_GME ||
gridInqType(gridID) == GRID_UNSTRUCTURED) )
if ( (operatorID == FILLMISS || operatorID == FILLMISSONESTEP) &&
(gridInqType(gridID) == GRID_GME || gridInqType(gridID) == GRID_UNSTRUCTURED) )
cdoAbort("%s data unsupported!", gridNamePtr(gridInqType(gridID)) );
field1.grid = gridID;
......@@ -495,9 +513,9 @@ void *Fillmiss(void *argument)
fill_method(&field1, &field2, nfill);
gridsize = gridInqSize(field2.grid);
nmiss = 0;
for ( i = 0; i < gridsize; ++i )
int gridsize = gridInqSize(field2.grid);
int nmiss = 0;
for ( int i = 0; i < gridsize; ++i )
if ( DBL_IS_EQUAL(field2.ptr[i], field2.missval) ) nmiss++;
streamWriteRecord(streamID2, field2.ptr, nmiss);
......
......@@ -408,7 +408,7 @@ void *Maggraph(void *argument);
#define SetgridOperators {"setgrid", "setgridtype", "setgridarea", "setgridmask", "unsetgridmask", "setgridnumber", "setgriduri"}
#define SethaloOperators {"sethalo", "tpnhalo"}
#define SetmissOperators {"setmissval", "setctomiss", "setmisstoc", "setrtomiss", "setvrange"}
#define SetmisstonnOperators {"setmisstonn"}
#define SetmisstonnOperators {"setmisstonn", "setmisstodis"}
#define Setpartab0Operators {"setpartab"}
#define SetpartabOperators {"setpartabc", "setpartabp", "setpartabn"}
#define SetrcanameOperators {"setrcaname"}
......
......@@ -65,32 +65,45 @@ void nbr_check_distance(unsigned num_neighbors, const int *restrict nbr_add, dou
}
}
static
double nbr_compute_weights(unsigned num_neighbors, const int *restrict src_grid_mask, int *restrict nbr_mask, const int *restrict nbr_add, double *restrict nbr_dist)
{
/* Compute weights based on inverse distance if mask is false, eliminate those points */
double dist_tot = 0.; /* sum of neighbor distances (for normalizing) */
for ( unsigned n = 0; n < num_neighbors; ++n )
if ( src_grid_mask )
{
for ( unsigned n = 0; n < num_neighbors; ++n )
{
nbr_mask[n] = FALSE;
if ( nbr_add[n] >= 0 )
if ( src_grid_mask[nbr_add[n]] )
{
nbr_dist[n] = ONE/nbr_dist[n];
dist_tot = dist_tot + nbr_dist[n];
nbr_mask[n] = TRUE;
}
}
}
else
{
// printf("tgt_cell_add %ld %ld %d %g\n", tgt_cell_add, n, nbr_add[n], nbr_dist[n]);
nbr_mask[n] = FALSE;
/* Uwe Schulzweida: check if nbr_add is valid */
if ( nbr_add[n] >= 0 )
if ( src_grid_mask[nbr_add[n]] )
{
nbr_dist[n] = ONE/nbr_dist[n];
dist_tot = dist_tot + nbr_dist[n];
nbr_mask[n] = TRUE;
}
for ( unsigned n = 0; n < num_neighbors; ++n )
{
nbr_mask[n] = FALSE;
if ( nbr_add[n] >= 0 )
{
nbr_dist[n] = ONE/nbr_dist[n];
dist_tot = dist_tot + nbr_dist[n];
nbr_mask[n] = TRUE;
}
}
}
return dist_tot;
}
static
unsigned nbr_normalize_weights(unsigned num_neighbors, double dist_tot, const int *restrict nbr_mask, int *restrict nbr_add, double *restrict nbr_dist)
{
/* Normalize weights and store the link */
......@@ -132,7 +145,7 @@ double get_search_radius(void)
#define MAX_SEARCH_CELLS 25
static
void grid_search_nbr_reg2d(struct gridsearch *gs, int num_neighbors, remapgrid_t *src_grid, int *restrict nbr_add, double *restrict nbr_dist,
double plat, double plon, const int *restrict src_grid_dims)
double plon, double plat, const int *restrict src_grid_dims)
{
/*
Output variables:
......@@ -279,8 +292,8 @@ void grid_search_nbr_reg2d(struct gridsearch *gs, int num_neighbors, remapgrid_t
}
} /* grid_search_nbr_reg2d */
static
void grid_search_nbr(struct gridsearch *gs, int num_neighbors, int *restrict nbr_add, double *restrict nbr_dist, double plat, double plon)
void grid_search_nbr(struct gridsearch *gs, int num_neighbors, int *restrict nbr_add, double *restrict nbr_dist, double plon, double plat)
{
/*
Output variables:
......@@ -435,16 +448,14 @@ void remap_distwgt_weights(unsigned num_neighbors, remapgrid_t *src_grid, remapg
/* Find nearest grid points on source grid and distances to each point */
if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
grid_search_nbr_reg2d(gs, num_neighbors, src_grid, nbr_add, nbr_dist,
plat, plon, src_grid->dims);
plon, plat, src_grid->dims);
else
grid_search_nbr(gs, num_neighbors, nbr_add, nbr_dist, plat, plon);
grid_search_nbr(gs, num_neighbors, nbr_add, nbr_dist, plon, plat);
/* Compute weights based on inverse distance if mask is false, eliminate those points */
double dist_tot = nbr_compute_weights(num_neighbors, src_grid->mask, nbr_mask, nbr_add, nbr_dist);
/* Normalize weights and store the link */
unsigned nadds = nbr_normalize_weights(num_neighbors, dist_tot, nbr_mask, nbr_add, nbr_dist);
for ( unsigned n = 0; n < nadds; ++n )
......@@ -541,16 +552,14 @@ void remap_distwgt(unsigned num_neighbors, remapgrid_t *src_grid, remapgrid_t *t
/* Find nearest grid points on source grid and distances to each point */
if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
grid_search_nbr_reg2d(gs, num_neighbors, src_grid, nbr_add, nbr_dist,
plat, plon, src_grid->dims);
plon, plat, src_grid->dims);
else
grid_search_nbr(gs, num_neighbors, nbr_add, nbr_dist, plat, plon);
grid_search_nbr(gs, num_neighbors, nbr_add, nbr_dist, plon, plat);
/* Compute weights based on inverse distance if mask is false, eliminate those points */
double dist_tot = nbr_compute_weights(num_neighbors, src_grid->mask, nbr_mask, nbr_add, nbr_dist);
/* Normalize weights and store the link */
unsigned nadds = nbr_normalize_weights(num_neighbors, dist_tot, nbr_mask, nbr_add, nbr_dist);
for ( unsigned n = 0; n < nadds; ++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