Commit 1c68b574 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

expr_var_con: added missval support

parent 75ee4e9f
......@@ -136,8 +136,8 @@ void *Expr(void *argument)
parse_arg.init = 1;
parse_arg.vlistID1 = vlistID1;
parse_arg.vlistID2 = vlistID2;
parse_arg.nvars1 = 0;
parse_arg.debug = 0;
parse_arg.nvars1 = 0;
parse_arg.debug = 0;
parse_arg.gridID2 = -1;
parse_arg.zaxisID2 = -1;
parse_arg.timeID2 = -1;
......@@ -150,8 +150,7 @@ void *Expr(void *argument)
parse_arg.init = 0;
nvars2 = vlistNvars(vlistID2);
if ( nvars2 == 0 ) cdoAbort("No variable in output!");
if ( nvars2 == 0 ) cdoAbort("No output variable found!");
if ( cdoVerbose ) vlistPrint(vlistID2);
......@@ -175,7 +174,7 @@ void *Expr(void *argument)
{
gridID = vlistInqVarGrid(vlistID1, varID);
zaxisID = vlistInqVarZaxis(vlistID1, varID);
/* parse_arg.missval = vlistInqVarMissval(vlistID1, varID); */
/* parse_arg.missval[varID] = vlistInqVarMissval(vlistID1, varID); */
gridsize = gridInqSize(gridID);
nlevel = zaxisInqSize(zaxisID);
......@@ -205,6 +204,8 @@ void *Expr(void *argument)
streamDefTimestep(streamID2, tsID);
for ( varID = 0; varID < nvars; varID++ ) parse_arg.nmiss[varID] = 0;
for ( recID = 0; recID < nrecs; recID++ )
{
streamInqRecord(streamID1, &varID, &levelID);
......@@ -215,11 +216,14 @@ void *Expr(void *argument)
offset = gridsize*levelID;
single1 = parse_arg.vardata1[varID] + offset;
streamReadRecord(streamID1, single1, &nmiss);
parse_arg.nmiss[varID] += nmiss;
/*
if ( nmiss && lwarn )
{
cdoWarning("Missing values unsupported for this operator!");
lwarn = FALSE;
}
*/
}
}
......@@ -244,23 +248,16 @@ void *Expr(void *argument)
gridsize = gridInqSize(gridID);
nlevel = zaxisInqSize(zaxisID);
/* nmiss = parse_arg.nmiss; */
/* if ( nmiss ) fprintf(stdout, "out nmiss = %d\n", nmiss); */
for ( levelID = 0; levelID < nlevel; levelID++ )
{
long i;
offset = gridsize*levelID;
single2 = parse_arg.vardata2[varID] + offset;
nmiss = 0;
if ( missval < -1.e30 || missval > 1.e30 )
{
int i;
for ( i = 0; i < gridsize; i++ )
if ( single2[i] < -1.e30 || single2[i] > 1.e30 )
{
single2[i] = missval;
nmiss++;
}
}
for ( i = 0; i < gridsize; i++ )
if ( DBL_IS_EQUAL(single2[i], missval) ) nmiss++;
streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, single2, nmiss);
}
......
......@@ -5,10 +5,14 @@
#include "cdi.h"
#include "cdo.h"
#include "cdo_int.h"
#include "field.h"
#include "expr.h"
#include "expr_yacc.h"
#define UMISS(func, oper) fprintf(stderr, "Internal problem in %s: missing value support not implemented for operation '%c'!\n", func, oper);
typedef struct {
int type;
char *name; /* function name */
......@@ -73,12 +77,15 @@ nodeType *expr_con_var(int oper, nodeType *p1, nodeType *p2)
nodeType *p;
int ngp, i;
int nlev, k;
int nmiss;
int gridID, zaxisID;
double missval;
double missval1, missval2;
gridID = p2->gridID;
zaxisID = p2->zaxisID;
missval = p2->missval;
gridID = p2->gridID;
zaxisID = p2->zaxisID;
nmiss = p2->nmiss;
missval1 = p2->missval;
missval2 = p2->missval;
ngp = gridInqSize(gridID);
nlev = zaxisInqSize(zaxisID);
......@@ -90,41 +97,87 @@ nodeType *expr_con_var(int oper, nodeType *p1, nodeType *p2)
p->u.var.nm = strdupx("tmp");
p->gridID = gridID;
p->zaxisID = zaxisID;
p->missval = missval;
p->missval = missval1;
p->data = (double *) malloc(ngp*nlev*sizeof(double));
switch ( oper )
{
case '+':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->u.con.value + p2->data[i+k*ngp];
if ( nmiss > 0 )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = ADD(p1->u.con.value, p2->data[i+k*ngp]);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->u.con.value + p2->data[i+k*ngp];
}
break;
case '-':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->u.con.value - p2->data[i+k*ngp];
if ( nmiss > 0 )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = SUB(p1->u.con.value, p2->data[i+k*ngp]);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->u.con.value - p2->data[i+k*ngp];
}
break;
case '*':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->u.con.value * p2->data[i+k*ngp];
if ( nmiss > 0 )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = MUL(p1->u.con.value, p2->data[i+k*ngp]);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->u.con.value * p2->data[i+k*ngp];
}
break;
case '/':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->u.con.value / p2->data[i+k*ngp];
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = DIV(p1->u.con.value, p2->data[i+k*ngp]);
}
break;
case '^':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->u.con.value, p2->data[i+k*ngp]);
if ( nmiss > 0 )
{
UMISS(func, oper);
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->u.con.value, p2->data[i+k*ngp]);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->u.con.value, p2->data[i+k*ngp]);
}
break;
default:
cdoAbort("%s: operator %c unsupported!", func, oper);
}
nmiss = 0;
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
if ( DBL_IS_EQUAL(p->data[i+k*ngp], missval1) ) nmiss++;
p->nmiss = nmiss;
if ( p2->tmpvar ) free(p2->data);
return (p);
......@@ -137,12 +190,15 @@ nodeType *expr_var_con(int oper, nodeType *p1, nodeType *p2)
nodeType *p;
int ngp, i;
int nlev, k;
int nmiss;
int gridID, zaxisID;
double missval;
double missval1, missval2;
gridID = p1->gridID;
zaxisID = p1->zaxisID;
missval = p1->missval;
gridID = p1->gridID;
zaxisID = p1->zaxisID;
nmiss = p1->nmiss;
missval1 = p1->missval;
missval2 = p1->missval;
ngp = gridInqSize(gridID);
nlev = zaxisInqSize(zaxisID);
......@@ -154,41 +210,94 @@ nodeType *expr_var_con(int oper, nodeType *p1, nodeType *p2)
p->u.var.nm = strdupx("tmp");
p->gridID = gridID;
p->zaxisID = zaxisID;
p->missval = missval;
p->missval = missval1;
p->data = (double *) malloc(ngp*nlev*sizeof(double));
switch ( oper )
{
case '+':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+k*ngp] + p2->u.con.value;
if ( nmiss > 0 )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = ADD(p1->data[i+k*ngp], p2->u.con.value);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+k*ngp] + p2->u.con.value;
}
break;
case '-':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+k*ngp] - p2->u.con.value;
if ( nmiss > 0 )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = SUB(p1->data[i+k*ngp], p2->u.con.value);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+k*ngp] - p2->u.con.value;
}
break;
case '*':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+k*ngp] * p2->u.con.value;
if ( nmiss > 0 )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = MUL(p1->data[i+k*ngp], p2->u.con.value);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+k*ngp] * p2->u.con.value;
}
break;
case '/':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+k*ngp] / p2->u.con.value;
if ( nmiss > 0 || IS_EQUAL(p2->u.con.value, 0) )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = DIV(p1->data[i+k*ngp], p2->u.con.value);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+k*ngp] / p2->u.con.value;
}
break;
case '^':
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->data[i+k*ngp], p2->u.con.value);
if ( nmiss > 0 )
{
UMISS(func, oper);
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->data[i+k*ngp], p2->u.con.value);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->data[i+k*ngp], p2->u.con.value);
}
break;
default:
cdoAbort("%s: operator %c unsupported!", func, oper);
}
nmiss = 0;
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
if ( DBL_IS_EQUAL(p->data[i+k*ngp], missval1) ) nmiss++;
p->nmiss = nmiss;
if ( p1->tmpvar ) free(p1->data);
return (p);
......@@ -305,6 +414,9 @@ void ex_copy(nodeType *p2, nodeType *p1)
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p2->data[i+k*ngp] = p1->data[i+k*ngp];
p2->missval = p1->missval;
p2->nmiss = p1->nmiss;
}
static
......@@ -584,7 +696,7 @@ nodeType *ex(nodeType *p, parse_parm_t *parse_arg)
if ( varID == nvars )
{
cdoAbort("variable >%s< not found!", p->u.var.nm);
cdoAbort("Variable >%s< not found!", p->u.var.nm);
}
else
{
......@@ -620,7 +732,10 @@ nodeType *ex(nodeType *p, parse_parm_t *parse_arg)
p->missval = missval;
p->nmiss = 0;
if ( ! parse_arg->init )
p->data = parse_arg->vardata1[varID];
{
p->data = parse_arg->vardata1[varID];
p->nmiss = parse_arg->nmiss[varID];
}
p->tmpvar = 0;
rnode = p;
}
......@@ -689,8 +804,8 @@ nodeType *ex(nodeType *p, parse_parm_t *parse_arg)
p->gridID = parse_arg->gridID2;
p->zaxisID = parse_arg->zaxisID2;
p->missval = missval;
p->data = parse_arg->vardata2[varID];
p->tmpvar = 0;
p->data = parse_arg->vardata2[varID];
p->tmpvar = 0;
ex_copy(p, rnode);
......
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