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

New operator: setmisstonn - Set missing value to nearest neightbour

parent 27c7a0b2
Hide whitespace changes
Inline Side-by-side
ChangeLog
View file @ c163187d
......@@ -3,6 +3,10 @@
* using CDI library version 1.7.0
* Version 1.7.0 released
2015-06-27 Uwe Schulzweida
* New operator: setmisstonn - Set missing value to nearest neightbour
2015-06-26 Uwe Schulzweida
* remapdis: fixed seg fault on blizzard with large target grid (possibly compiler bug)
......
NEWS
View file @ c163187d
......@@ -4,6 +4,7 @@ CDO NEWS
Version 1.7.0 (28 October 2015):
New operators:
* setmisstonn: Set missing value to nearest neightbour
* vertstd1: Vertical standard deviation [Divisor is (n-1)]
* vertvar1: Vertical variance [Divisor is (n-1)]
* seasvar1: Seasonal variance [Divisor is (n-1)]
......
OPERATORS
View file @ c163187d
......@@ -163,6 +163,7 @@ Operator catalog:
Setmiss setmisstoc Set missing value to constant
Setmiss setrtomiss Set range to missing value
Setmiss setvrange Set valid range
Setmiss setmisstonn Set missing value to nearest neightbour
-------------------------------------------------------------
Arithmetic
-------------------------------------------------------------
......
doc/cdo_refcard.pdf
View file @ c163187d
No preview for this file type
doc/tex/mod/Setmiss
View file @ c163187d
......@@ -4,7 +4,7 @@
@Title = Set missing value
@Section = Modification
@Arguments = ifile ofile
@Operators = setmissval setctomiss setmisstoc setrtomiss setvrange
@Operators = setmissval setctomiss setmisstoc setrtomiss setvrange setmisstonn
@BeginDescription
This module sets part of a field to missing value or missing values
......@@ -129,6 +129,30 @@ o(t,x) = \left\{
@EndOperator
@BeginOperator_setmisstonn
@Title = Set missing value to nearest neightbour
@BeginDescription
Set all missing values to the nearest non missing value.
@IfMan
/ i(t,y) if i(t,x) EQ miss AND i(t,y) NE miss
o(t,x) =
\ i(t,x) if i(t,x) NE miss
@EndifMan
@IfDoc
@BeginMath
o(t,x) = \left\{
\begin{array}{cll}
i(t,y) & \mbox{if} \;\; i(t,x) \geq \mbox{miss} \wedge i(t,y) \neq \mbox{miss} \\
i(t,x) & \mbox{if} \;\; i(t,x) \neq \mbox{miss} \\
\end{array} \right.
@EndMath
@EndifDoc
@EndDescription
@EndOperator
@BeginParameter rmax
@Item = newmiss
FLOAT New missing value
......
src/Fillmiss.c
View file @ c163187d
......@@ -280,7 +280,7 @@ static double dist_sq( double *a1, double *a2, int dims ) {
}
void fillmiss_nn(field_t *field1, field_t *field2, int maxfill)
void setmisstonn(field_t *field1, field_t *field2, int maxfill)
{
int gridID = field1->grid;
double missval = field1->missval;
......@@ -359,7 +359,7 @@ void fillmiss_nn(field_t *field1, field_t *field2, int maxfill)
finish = clock();
printf("kd_tree created: %.2f seconds\n", ((double)(finish-start))/CLOCKS_PER_SEC);
if ( cdoVerbose ) printf("kd_tree created: %.2f seconds\n", ((double)(finish-start))/CLOCKS_PER_SEC);
start = clock();
......@@ -409,7 +409,7 @@ void fillmiss_nn(field_t *field1, field_t *field2, int maxfill)
*/
finish = clock();
printf("kd_tree nearest: %.2f seconds\n", ((double)(finish-start))/CLOCKS_PER_SEC);
if ( cdoVerbose ) printf("kd_tree nearest: %.2f seconds\n", ((double)(finish-start))/CLOCKS_PER_SEC);
kd_free(pointTree);
......@@ -430,7 +430,7 @@ void *Fillmiss(void *argument)
int FILLMISS = cdoOperatorAdd("fillmiss" , 0, 0, "nfill");
int FILLMISSONESTEP = cdoOperatorAdd("fillmiss2" , 0, 0, "nfill");
int FILLMISSNN = cdoOperatorAdd("fillmissnn" , 0, 0, "nfill");
int SETMISSTONN = cdoOperatorAdd("setmisstonn" , 0, 0, "nfill");
int operatorID = cdoOperatorID();
......@@ -442,9 +442,9 @@ void *Fillmiss(void *argument)
{
fill_method = &fillmiss_one_step;
}
else if ( operatorID == FILLMISSNN )
else if ( operatorID == SETMISSTONN )
{
fill_method = &fillmiss_nn;
fill_method = &setmisstonn;
}
/* Argument handling */
......
src/kdtree.c
View file @ c163187d
......@@ -25,6 +25,7 @@ IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
OF SUCH DAMAGE.
*/
/* single nearest neighbor search written by Tamas Nepusz <tamas@cs.rhul.ac.uk> */
/* U. Schulzweida, 20150627: removed float interface, set MAX_DIM=3 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
......@@ -35,28 +36,16 @@ OF SUCH DAMAGE.
#include <malloc.h>
#endif
#define NO_ALLOCA
#define MAX_DIM 3
#ifdef USE_LIST_NODE_ALLOCATOR
#ifndef NO_PTHREADS
#include <pthread.h>
#else
#ifndef I_WANT_THREAD_BUGS
#error "You are compiling with the fast list node allocator, with pthreads disabled! This WILL break if used from multiple threads."
#endif /* I want thread bugs */
#endif /* pthread support */
#endif /* use list node allocator */
struct kdhyperrect {
int dim;
double *min, *max; /* minimum/maximum coords */
double min[MAX_DIM], max[MAX_DIM]; /* minimum/maximum coords */
};
struct kdnode {
double *pos;
double pos[MAX_DIM];
int dir;
void *data;
......@@ -96,14 +85,8 @@ static struct kdhyperrect* hyperrect_duplicate(const struct kdhyperrect *rect);
static void hyperrect_extend(struct kdhyperrect *rect, const double *pos);
static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos);
#ifdef USE_LIST_NODE_ALLOCATOR
static struct res_node *alloc_resnode(void);
static void free_resnode(struct res_node*);
#else
#define alloc_resnode() malloc(sizeof(struct res_node))
#define free_resnode(n) free(n)
#endif
struct kdtree *kd_create(int k)
......@@ -140,7 +123,6 @@ static void clear_rec(struct kdnode *node, void (*destr)(void*))
if(destr) {
destr(node->data);
}
free(node->pos);
free(node);
}
......@@ -170,10 +152,6 @@ static int insert_rec(struct kdnode **nptr, const double *pos, void *data, int d
if(!(node = malloc(sizeof *node))) {
return -1;
}
if(!(node->pos = malloc(dim * sizeof *node->pos))) {
free(node);
return -1;
}
memcpy(node->pos, pos, dim * sizeof *node->pos);
node->data = data;
node->dir = dir;
......@@ -205,38 +183,6 @@ int kd_insert(struct kdtree *tree, const double *pos, void *data)
return 0;
}
int kd_insertf(struct kdtree *tree, const float *pos, void *data)
{
static double sbuf[16];
double *bptr, *buf = 0;
int res, dim = tree->dim;
if(dim > 16) {
#ifndef NO_ALLOCA
if(dim <= 256)
bptr = buf = alloca(dim * sizeof *bptr);
else
#endif
if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
return -1;
}
} else {
bptr = buf = sbuf;
}
while(dim-- > 0) {
*bptr++ = *pos++;
}
res = kd_insert(tree, buf, data);
#ifndef NO_ALLOCA
if(tree->dim > 256)
#else
if(tree->dim > 16)
#endif
free(buf);
return res;
}
int kd_insert3(struct kdtree *tree, double x, double y, double z, void *data)
{
......@@ -247,16 +193,8 @@ int kd_insert3(struct kdtree *tree, double x, double y, double z, void *data)
return kd_insert(tree, buf, data);
}
int kd_insert3f(struct kdtree *tree, float x, float y, float z, void *data)
{
double buf[3];
buf[0] = x;
buf[1] = y;
buf[2] = z;
return kd_insert(tree, buf, data);
}
static int find_nearest(struct kdnode *node, const double *pos, double range, struct res_node *list, int ordered, int dim)
static
int find_nearest(struct kdnode *node, const double *pos, double range, struct res_node *list, int ordered, int dim)
{
double dist_sq, dx;
int i, ret, added_res = 0;
......@@ -453,40 +391,6 @@ struct kdres *kd_nearest(struct kdtree *kd, const double *pos)
}
}
struct kdres *kd_nearestf(struct kdtree *tree, const float *pos)
{
static double sbuf[16];
double *bptr, *buf = 0;
int dim = tree->dim;
struct kdres *res;
if(dim > 16) {
#ifndef NO_ALLOCA
if(dim <= 256)
bptr = buf = alloca(dim * sizeof *bptr);
else
#endif
if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
return 0;
}
} else {
bptr = buf = sbuf;
}
while(dim-- > 0) {
*bptr++ = *pos++;
}
res = kd_nearest(tree, buf);
#ifndef NO_ALLOCA
if(tree->dim > 256)
#else
if(tree->dim > 16)
#endif
free(buf);
return res;
}
struct kdres *kd_nearest3(struct kdtree *tree, double x, double y, double z)
{
double pos[3];
......@@ -496,15 +400,6 @@ struct kdres *kd_nearest3(struct kdtree *tree, double x, double y, double z)
return kd_nearest(tree, pos);
}
struct kdres *kd_nearest3f(struct kdtree *tree, float x, float y, float z)
{
double pos[3];
pos[0] = x;
pos[1] = y;
pos[2] = z;
return kd_nearest(tree, pos);
}
/* ---- nearest N search ---- */
/*
static kdres *kd_nearest_n(struct kdtree *kd, const double *pos, int num)
......@@ -555,40 +450,6 @@ struct kdres *kd_nearest_range(struct kdtree *kd, const double *pos, double rang
return rset;
}
struct kdres *kd_nearest_rangef(struct kdtree *kd, const float *pos, float range)
{
static double sbuf[16];
double *bptr, *buf = 0;
int dim = kd->dim;
struct kdres *res;
if(dim > 16) {
#ifndef NO_ALLOCA
if(dim <= 256)
bptr = buf = alloca(dim * sizeof *bptr);
else
#endif
if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
return 0;
}
} else {
bptr = buf = sbuf;
}
while(dim-- > 0) {
*bptr++ = *pos++;
}
res = kd_nearest_range(kd, buf, range);
#ifndef NO_ALLOCA
if(kd->dim > 256)
#else
if(kd->dim > 16)
#endif
free(buf);
return res;
}
struct kdres *kd_nearest_range3(struct kdtree *tree, double x, double y, double z, double range)
{
double buf[3];
......@@ -598,15 +459,6 @@ struct kdres *kd_nearest_range3(struct kdtree *tree, double x, double y, double
return kd_nearest_range(tree, buf, range);
}
struct kdres *kd_nearest_range3f(struct kdtree *tree, float x, float y, float z, float range)
{
double buf[3];
buf[0] = x;
buf[1] = y;
buf[2] = z;
return kd_nearest_range(tree, buf, range);
}
void kd_res_free(struct kdres *rset)
{
clear_results(rset);
......@@ -646,20 +498,6 @@ void *kd_res_item(struct kdres *rset, double *pos)
return 0;
}
void *kd_res_itemf(struct kdres *rset, float *pos)
{
if(rset->riter) {
if(pos) {
int i;
for(i=0; i<rset->tree->dim; i++) {
pos[i] = rset->riter->item->pos[i];
}
}
return rset->riter->item->data;
}
return 0;
}
void *kd_res_item3(struct kdres *rset, double *x, double *y, double *z)
{
if(rset->riter) {
......@@ -670,16 +508,6 @@ void *kd_res_item3(struct kdres *rset, double *x, double *y, double *z)
return 0;
}
void *kd_res_item3f(struct kdres *rset, float *x, float *y, float *z)
{
if(rset->riter) {
if(*x) *x = rset->riter->item->pos[0];
if(*y) *y = rset->riter->item->pos[1];
if(*z) *z = rset->riter->item->pos[2];
}
return 0;
}
void *kd_res_item_data(struct kdres *set)
{
return kd_res_item(set, 0);
......@@ -696,15 +524,6 @@ static struct kdhyperrect* hyperrect_create(int dim, const double *min, const do
}
rect->dim = dim;
if (!(rect->min = malloc(size))) {
free(rect);
return 0;
}
if (!(rect->max = malloc(size))) {
free(rect->min);
free(rect);
return 0;
}
memcpy(rect->min, min, size);
memcpy(rect->max, max, size);
......@@ -713,8 +532,6 @@ static struct kdhyperrect* hyperrect_create(int dim, const double *min, const do
static void hyperrect_free(struct kdhyperrect *rect)
{
free(rect->min);
free(rect->max);
free(rect);
}
......@@ -753,54 +570,6 @@ static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos)
return result;
}
/* ---- static helpers ---- */
#ifdef USE_LIST_NODE_ALLOCATOR
/* special list node allocators. */
static struct res_node *free_nodes;
#ifndef NO_PTHREADS
static pthread_mutex_t alloc_mutex = PTHREAD_MUTEX_INITIALIZER;
#endif
static struct res_node *alloc_resnode(void)
{
struct res_node *node;
#ifndef NO_PTHREADS
pthread_mutex_lock(&alloc_mutex);
#endif
if(!free_nodes) {
node = malloc(sizeof *node);
} else {
node = free_nodes;
free_nodes = free_nodes->next;
node->next = 0;
}
#ifndef NO_PTHREADS
pthread_mutex_unlock(&alloc_mutex);
#endif
return node;
}
static void free_resnode(struct res_node *node)
{
#ifndef NO_PTHREADS
pthread_mutex_lock(&alloc_mutex);
#endif
node->next = free_nodes;
free_nodes = node;
#ifndef NO_PTHREADS
pthread_mutex_unlock(&alloc_mutex);
#endif
}
#endif /* list node allocator or not */
/* inserts the item. if dist_sq is >= 0, then do an ordered insert */
/* TODO make the ordering code use heapsort */
......
src/kdtree.h
View file @ c163187d
......@@ -52,18 +52,14 @@ void kd_data_destructor(struct kdtree *tree, void (*destr)(void*));
/* insert a node, specifying its position, and optional data */
int kd_insert(struct kdtree *tree, const double *pos, void *data);
int kd_insertf(struct kdtree *tree, const float *pos, void *data);
int kd_insert3(struct kdtree *tree, double x, double y, double z, void *data);
int kd_insert3f(struct kdtree *tree, float x, float y, float z, void *data);
/* Find the nearest node from a given point.
*
* This function returns a pointer to a result set with at most one element.
*/
struct kdres *kd_nearest(struct kdtree *tree, const double *pos);
struct kdres *kd_nearestf(struct kdtree *tree, const float *pos);
struct kdres *kd_nearest3(struct kdtree *tree, double x, double y, double z);
struct kdres *kd_nearest3f(struct kdtree *tree, float x, float y, float z);
/* Find the N nearest nodes from a given point.
*
......@@ -75,9 +71,7 @@ struct kdres *kd_nearest3f(struct kdtree *tree, float x, float y, float z);
*/
/*
struct kdres *kd_nearest_n(struct kdtree *tree, const double *pos, int num);
struct kdres *kd_nearest_nf(struct kdtree *tree, const float *pos, int num);
struct kdres *kd_nearest_n3(struct kdtree *tree, double x, double y, double z);
struct kdres *kd_nearest_n3f(struct kdtree *tree, float x, float y, float z);
*/
/* Find any nearest nodes from a given point within a range.
......@@ -89,9 +83,7 @@ struct kdres *kd_nearest_n3f(struct kdtree *tree, float x, float y, float z);
* The result set must be deallocated with kd_res_free after use.
*/
struct kdres *kd_nearest_range(struct kdtree *tree, const double *pos, double range);
struct kdres *kd_nearest_rangef(struct kdtree *tree, const float *pos, float range);
struct kdres *kd_nearest_range3(struct kdtree *tree, double x, double y, double z, double range);
struct kdres *kd_nearest_range3f(struct kdtree *tree, float x, float y, float z, float range);
/* frees a result set returned by kd_nearest_range() */
void kd_res_free(struct kdres *set);
......@@ -114,9 +106,7 @@ int kd_res_next(struct kdres *set);
* and optionally sets its position to the pointers(s) if not null.
*/
void *kd_res_item(struct kdres *set, double *pos);
void *kd_res_itemf(struct kdres *set, float *pos);
void *kd_res_item3(struct kdres *set, double *x, double *y, double *z);
void *kd_res_item3f(struct kdres *set, float *x, float *y, float *z);
/* equivalent to kd_res_item(set, 0) */
void *kd_res_item_data(struct kdres *set);
......
src/modules.c
View file @ c163187d
......@@ -318,7 +318,7 @@ void *Maggraph(void *argument);
#define FCOperators {"fc2sp", "sp2fc", "fc2gp", "gp2fc"}
#define FiledesOperators {"filedes", "griddes", "griddes2", "zaxisdes", "vct", "vct2", "pardes", \
"vlist", "partab", "partab2", "spartab"}
#define FillmissOperators {"fillmiss", "fillmiss2", "fillmissnn"}
#define FillmissOperators {"fillmiss", "fillmiss2"}
#define FilterOperators {"bandpass", "highpass", "lowpass"}
#define FldrmsOperators {"fldrms"}
#define FldstatOperators {"fldmin", "fldmax", "fldsum", "fldmean", "fldavg", "fldstd", "fldstd1", "fldvar", "fldvar1", "fldpctl"}
......@@ -400,6 +400,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 Setpartab0Operators {"setpartab"}
#define SetpartabOperators {"setpartabc", "setpartabp", "setpartabn"}
#define SetrcanameOperators {"setrcaname"}
......@@ -665,6 +666,7 @@ static modules_t Modules[] =
{ Setgrid, SetgridHelp, SetgridOperators, CDI_BOTH, 1, 1 },
{ Sethalo, SethaloHelp, SethaloOperators, CDI_REAL, 1, 1 },
{ Setmiss, SetmissHelp, SetmissOperators, CDI_REAL, 1, 1 },
{ Fillmiss, SetmissHelp, SetmisstonnOperators, CDI_REAL, 1, 1 },
{ Setpartab, SetHelp, Setpartab0Operators, CDI_REAL, 1, 1 },
{ Setpartab, SetpartabHelp, SetpartabOperators, CDI_REAL, 1, 1 },
{ Setrcaname, NULL, SetrcanameOperators, CDI_REAL, 1, 1 },
......
src/operator_help.h
View file @ c163187d
......@@ -1226,7 +1226,8 @@ static char *EnlargeHelp[] = {
static char *SetmissHelp[] = {
"NAME",
" setmissval, setctomiss, setmisstoc, setrtomiss, setvrange - Set missing value",
" setmissval, setctomiss, setmisstoc, setrtomiss, setvrange, setmisstonn - ",
" Set missing value",
"",
"SYNOPSIS",
" setmissval,newmiss ifile ofile",
......@@ -1234,6 +1235,7 @@ static char *SetmissHelp[] = {
" setmisstoc,c ifile ofile",
" setrtomiss,rmin,rmax ifile ofile",
" setvrange,rmin,rmax ifile ofile",
" setmisstonn ifile ofile",
"",
"DESCRIPTION",
" This module sets part of a field to missing value or missing values",
......@@ -1241,26 +1243,31 @@ static char *SetmissHelp[] = {
" chosen operator.",
"",
"OPERATORS",
" setmissval Set a new missing value",
" / newmiss if i(t,x) EQ miss",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) NE miss",
" setctomiss Set constant to missing value",
" / miss if i(t,x) EQ c",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) NE c",
" setmisstoc Set missing value to constant",
" / c if i(t,x) EQ miss",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) NE miss",
" setrtomiss Set range to missing value",
" / miss if i(t,x) GE rmin AND i(t,x) LE rmax",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) LT rmin OR i(t,x) GT rmax",
" setvrange Set valid range",
" / miss if i(t,x) LT rmin OR i(t,x) GT rmax",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) GE rmin AND i(t,x) LE rmax",
" setmissval Set a new missing value",
" / newmiss if i(t,x) EQ miss",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) NE miss",
" setctomiss Set constant to missing value",
" / miss if i(t,x) EQ c",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) NE c",
" setmisstoc Set missing value to constant",
" / c if i(t,x) EQ miss",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) NE miss",
" setrtomiss Set range to missing value",
" / miss if i(t,x) GE rmin AND i(t,x) LE rmax",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) LT rmin OR i(t,x) GT rmax",
" setvrange Set valid range",
" / miss if i(t,x) LT rmin OR i(t,x) GT rmax",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) GE rmin AND i(t,x) LE rmax",
" setmisstonn Set missing value to nearest neightbour",
" Set all missing values to the nearest non missing value.",
" / i(t,y) if i(t,x) EQ miss AND i(t,y) NE miss",
" o(t,x) = ",
" \\ i(t,x) if i(t,x) NE miss",
"",
"PARAMETER",
" newmiss FLOAT New missing value",
......
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