Commit 4aecbc8f authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added new operator timselrange.

parent c3c4e7ea
...@@ -6,6 +6,7 @@ ...@@ -6,6 +6,7 @@
2017-05-30 Uwe Schulzweida 2017-05-30 Uwe Schulzweida
* New operator seasrange: seasonal range (seasmax-seasmin) * New operator seasrange: seasonal range (seasmax-seasmin)
* New operator timselrange: time selection range (timmax-timmin)
2017-05-26 Uwe Schulzweida 2017-05-26 Uwe Schulzweida
......
...@@ -159,20 +159,20 @@ void *Seasstat(void *argument) ...@@ -159,20 +159,20 @@ void *Seasstat(void *argument)
recinfo[recID].levelID = levelID; recinfo[recID].levelID = levelID;
} }
field_type *pvar1 = &vars1[varID][levelID]; field_type *pvars1 = &vars1[varID][levelID];
field_type *pvars2 = vars2 ? &vars2[varID][levelID] : NULL;
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID)); gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
if ( nsets == 0 ) if ( nsets == 0 )
{ {
streamReadRecord(streamID1, pvar1->ptr, &nmiss); streamReadRecord(streamID1, pvars1->ptr, &nmiss);
pvar1->nmiss = (size_t)nmiss; pvars1->nmiss = (size_t)nmiss;
if ( lrange ) if ( lrange )
{ {
field_type *pvar2 = &vars2[varID][levelID]; pvars2->nmiss = (size_t)nmiss;
for ( int i = 0; i < gridsize; i++ ) for ( int i = 0; i < gridsize; i++ )
pvar2->ptr[i] = pvar1->ptr[i]; pvars2->ptr[i] = pvars1->ptr[i];
pvar2->nmiss = (size_t)nmiss;
} }
if ( nmiss > 0 || samp1[varID][levelID].ptr ) if ( nmiss > 0 || samp1[varID][levelID].ptr )
...@@ -181,15 +181,15 @@ void *Seasstat(void *argument) ...@@ -181,15 +181,15 @@ void *Seasstat(void *argument)
samp1[varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double)); samp1[varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
for ( int i = 0; i < gridsize; i++ ) for ( int i = 0; i < gridsize; i++ )
samp1[varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvar1->ptr[i], pvar1->missval); samp1[varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
} }
} }
else else
{ {
streamReadRecord(streamID1, field.ptr, &nmiss); streamReadRecord(streamID1, field.ptr, &nmiss);
field.nmiss = (size_t)nmiss; field.nmiss = (size_t)nmiss;
field.grid = pvar1->grid; field.grid = pvars1->grid;
field.missval = pvar1->missval; field.missval = pvars1->missval;
if ( field.nmiss > 0 || samp1[varID][levelID].ptr ) if ( field.nmiss > 0 || samp1[varID][levelID].ptr )
{ {
...@@ -201,25 +201,23 @@ void *Seasstat(void *argument) ...@@ -201,25 +201,23 @@ void *Seasstat(void *argument)
} }
for ( int i = 0; i < gridsize; i++ ) for ( int i = 0; i < gridsize; i++ )
if ( !DBL_IS_EQUAL(field.ptr[i], pvar1->missval) ) if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
samp1[varID][levelID].ptr[i]++; samp1[varID][levelID].ptr[i]++;
} }
if ( lvarstd ) if ( lvarstd )
{ {
field_type *pvar2 = &vars2[varID][levelID]; farsumq(pvars2, field);
farsumq(pvar2, field); farsum(pvars1, field);
farsum(pvar1, field);
} }
else if ( lrange ) else if ( lrange )
{ {
field_type *pvar2 = &vars2[varID][levelID]; farmin(pvars2, field);
farmin(pvar2, field); farmax(pvars1, field);
farmax(pvar1, field);
} }
else else
{ {
farfun(pvar1, field, operfunc); farfun(pvars1, field, operfunc);
} }
} }
} }
...@@ -229,12 +227,12 @@ void *Seasstat(void *argument) ...@@ -229,12 +227,12 @@ void *Seasstat(void *argument)
{ {
int varID = recinfo[recID].varID; int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID; int levelID = recinfo[recID].levelID;
field_type *pvar1 = &vars1[varID][levelID]; field_type *pvars1 = &vars1[varID][levelID];
field_type *pvar2 = &vars2[varID][levelID]; field_type *pvars2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue; if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
farmoq(pvar2, *pvar1); farmoq(pvars2, *pvars1);
} }
vdate1 = vdate; vdate1 = vdate;
...@@ -250,34 +248,34 @@ void *Seasstat(void *argument) ...@@ -250,34 +248,34 @@ void *Seasstat(void *argument)
{ {
int varID = recinfo[recID].varID; int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID; int levelID = recinfo[recID].levelID;
field_type *pvar1 = &vars1[varID][levelID]; field_type *pvars1 = &vars1[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue; if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
if ( samp1[varID][levelID].ptr == NULL ) if ( samp1[varID][levelID].ptr == NULL )
farcdiv(pvar1, (double)nsets); farcdiv(pvars1, (double)nsets);
else else
fardiv(pvar1, samp1[varID][levelID]); fardiv(pvars1, samp1[varID][levelID]);
} }
else if ( lvarstd ) else if ( lvarstd )
for ( int recID = 0; recID < maxrecs; recID++ ) for ( int recID = 0; recID < maxrecs; recID++ )
{ {
int varID = recinfo[recID].varID; int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID; int levelID = recinfo[recID].levelID;
field_type *pvar1 = &vars1[varID][levelID]; field_type *pvars1 = &vars1[varID][levelID];
field_type *pvar2 = &vars2[varID][levelID]; field_type *pvars2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue; if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
if ( samp1[varID][levelID].ptr == NULL ) if ( samp1[varID][levelID].ptr == NULL )
{ {
if ( lstd ) farcstd(pvar1, *pvar2, nsets, divisor); if ( lstd ) farcstd(pvars1, *pvars2, nsets, divisor);
else farcvar(pvar1, *pvar2, nsets, divisor); else farcvar(pvars1, *pvars2, nsets, divisor);
} }
else else
{ {
if ( lstd ) farstd(pvar1, *pvar2, samp1[varID][levelID], divisor); if ( lstd ) farstd(pvars1, *pvars2, samp1[varID][levelID], divisor);
else farvar(pvar1, *pvar2, samp1[varID][levelID], divisor); else farvar(pvars1, *pvars2, samp1[varID][levelID], divisor);
} }
} }
else if ( lrange ) else if ( lrange )
...@@ -285,12 +283,12 @@ void *Seasstat(void *argument) ...@@ -285,12 +283,12 @@ void *Seasstat(void *argument)
{ {
int varID = recinfo[recID].varID; int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID; int levelID = recinfo[recID].levelID;
field_type *pvar1 = &vars1[varID][levelID]; field_type *pvars1 = &vars1[varID][levelID];
field_type *pvar2 = &vars2[varID][levelID]; field_type *pvars2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue; if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
farsub(pvar1, *pvar2); farsub(pvars1, *pvars2);
} }
if ( cdoVerbose ) if ( cdoVerbose )
...@@ -320,12 +318,12 @@ void *Seasstat(void *argument) ...@@ -320,12 +318,12 @@ void *Seasstat(void *argument)
{ {
int varID = recinfo[recID].varID; int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID; int levelID = recinfo[recID].levelID;
field_type *pvar1 = &vars1[varID][levelID]; field_type *pvars1 = &vars1[varID][levelID];
if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue; if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
streamDefRecord(streamID2, varID, levelID); streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, pvar1->ptr, (int)pvar1->nmiss); streamWriteRecord(streamID2, pvars1->ptr, (int)pvars1->nmiss);
} }
if ( nrecs == 0 ) break; if ( nrecs == 0 ) break;
...@@ -337,11 +335,11 @@ void *Seasstat(void *argument) ...@@ -337,11 +335,11 @@ void *Seasstat(void *argument)
field_free(samp1, vlistID1); field_free(samp1, vlistID1);
if ( lvarstd ) field_free(vars2, vlistID1); if ( lvarstd ) field_free(vars2, vlistID1);
if ( field.ptr ) Free(field.ptr); dtlist_delete(dtlist);
Free(recinfo); Free(recinfo);
dtlist_delete(dtlist); if ( field.ptr ) Free(field.ptr);
streamClose(streamID2); streamClose(streamID2);
streamClose(streamID1); streamClose(streamID1);
......
...@@ -69,35 +69,35 @@ void *Timcumsum(void *argument) ...@@ -69,35 +69,35 @@ void *Timcumsum(void *argument)
{ {
streamInqRecord(streamID1, &varID, &levelID); streamInqRecord(streamID1, &varID, &levelID);
field_type *pvar1 = &vars1[varID][levelID]; field_type *pvars1 = &vars1[varID][levelID];
gridsize = gridInqSize(pvar1->grid); gridsize = gridInqSize(pvars1->grid);
if ( tsID == 0 ) if ( tsID == 0 )
{ {
streamReadRecord(streamID1, pvar1->ptr, &nmiss); streamReadRecord(streamID1, pvars1->ptr, &nmiss);
// pvar1->nmiss = (size_t)nmiss; // pvars1->nmiss = (size_t)nmiss;
if ( nmiss ) if ( nmiss )
for ( int i = 0; i < gridsize; ++i ) for ( int i = 0; i < gridsize; ++i )
if ( DBL_IS_EQUAL(pvar1->ptr[i], pvar1->missval) ) pvar1->ptr[i] = 0; if ( DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval) ) pvars1->ptr[i] = 0;
} }
else else
{ {
streamReadRecord(streamID1, field.ptr, &nmiss); streamReadRecord(streamID1, field.ptr, &nmiss);
// field.nmiss = (size_t)nmiss; // field.nmiss = (size_t)nmiss;
field.size = gridsize; field.size = gridsize;
field.grid = pvar1->grid; field.grid = pvars1->grid;
field.missval = pvar1->missval; field.missval = pvars1->missval;
if ( nmiss ) if ( nmiss )
for ( int i = 0; i < gridsize; ++i ) for ( int i = 0; i < gridsize; ++i )
if ( DBL_IS_EQUAL(field.ptr[i], pvar1->missval) ) field.ptr[i] = 0; if ( DBL_IS_EQUAL(field.ptr[i], pvars1->missval) ) field.ptr[i] = 0;
farfun(pvar1, field, func_sum); farfun(pvars1, field, func_sum);
} }
streamDefRecord(streamID2, varID, levelID); streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, pvar1->ptr, (int)pvar1->nmiss); streamWriteRecord(streamID2, pvars1->ptr, (int)pvars1->nmiss);
} }
tsID++; tsID++;
......
...@@ -18,6 +18,7 @@ ...@@ -18,6 +18,7 @@
/* /*
This module contains the following operators: This module contains the following operators:
Timselstat timselrange Time range range
Timselstat timselmin Time range minimum Timselstat timselmin Time range minimum
Timselstat timselmax Time range maximum Timselstat timselmax Time range maximum
Timselstat timselsum Time range sum Timselstat timselsum Time range sum
...@@ -36,6 +37,12 @@ ...@@ -36,6 +37,12 @@
#include "pstream.h" #include "pstream.h"
typedef struct {
short varID;
short levelID;
} recinfo_t;
void *Timselstat(void *argument) void *Timselstat(void *argument)
{ {
int timestat_date = TIMESTAT_MEAN; int timestat_date = TIMESTAT_MEAN;
...@@ -44,19 +51,19 @@ void *Timselstat(void *argument) ...@@ -44,19 +51,19 @@ void *Timselstat(void *argument)
int tsID; int tsID;
int nsets; int nsets;
int nmiss; int nmiss;
int nlevel;
cdoInitialize(argument); cdoInitialize(argument);
cdoOperatorAdd("timselmin", func_min, 0, NULL); cdoOperatorAdd("timselrange", func_range, 0, NULL);
cdoOperatorAdd("timselmax", func_max, 0, NULL); cdoOperatorAdd("timselmin", func_min, 0, NULL);
cdoOperatorAdd("timselsum", func_sum, 0, NULL); cdoOperatorAdd("timselmax", func_max, 0, NULL);
cdoOperatorAdd("timselmean", func_mean, 0, NULL); cdoOperatorAdd("timselsum", func_sum, 0, NULL);
cdoOperatorAdd("timselavg", func_avg, 0, NULL); cdoOperatorAdd("timselmean", func_mean, 0, NULL);
cdoOperatorAdd("timselvar", func_var, 0, NULL); cdoOperatorAdd("timselavg", func_avg, 0, NULL);
cdoOperatorAdd("timselvar1", func_var1, 0, NULL); cdoOperatorAdd("timselvar", func_var, 0, NULL);
cdoOperatorAdd("timselstd", func_std, 0, NULL); cdoOperatorAdd("timselvar1", func_var1, 0, NULL);
cdoOperatorAdd("timselstd1", func_std1, 0, NULL); cdoOperatorAdd("timselstd", func_std, 0, NULL);
cdoOperatorAdd("timselstd1", func_std1, 0, NULL);
int operatorID = cdoOperatorID(); int operatorID = cdoOperatorID();
int operfunc = cdoOperatorF1(operatorID); int operfunc = cdoOperatorF1(operatorID);
...@@ -65,16 +72,16 @@ void *Timselstat(void *argument) ...@@ -65,16 +72,16 @@ void *Timselstat(void *argument)
int nargc = operatorArgc(); int nargc = operatorArgc();
int ndates = parameter2int(operatorArgv()[0]); int ndates = parameter2int(operatorArgv()[0]);
int noffset = 0, nskip = 0; int noffset = (nargc > 1) ? parameter2int(operatorArgv()[1]) : 0;
if ( nargc > 1 ) noffset = parameter2int(operatorArgv()[1]); int nskip = (nargc > 2) ? parameter2int(operatorArgv()[2]) : 0;
if ( nargc > 2 ) nskip = parameter2int(operatorArgv()[2]);
if ( cdoVerbose ) cdoPrint("nsets = %d, noffset = %d, nskip = %d", ndates, noffset, nskip); if ( cdoVerbose ) cdoPrint("nsets = %d, noffset = %d, nskip = %d", ndates, noffset, nskip);
int lmean = operfunc == func_mean || operfunc == func_avg; bool lrange = operfunc == func_range;
int lstd = operfunc == func_std || operfunc == func_std1; bool lmean = operfunc == func_mean || operfunc == func_avg;
int lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1; bool lstd = operfunc == func_std || operfunc == func_std1;
int divisor = operfunc == func_std1 || operfunc == func_var1; bool lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
int streamID1 = streamOpenRead(cdoStreamName(0)); int streamID1 = streamOpenRead(cdoStreamName(0));
...@@ -89,11 +96,9 @@ void *Timselstat(void *argument) ...@@ -89,11 +96,9 @@ void *Timselstat(void *argument)
streamDefVlist(streamID2, vlistID2); streamDefVlist(streamID2, vlistID2);
int nvars = vlistNvars(vlistID1); int maxrecs = vlistNrecs(vlistID1);
int nrecords = vlistNrecs(vlistID1);
int *recVarID = (int*) Malloc(nrecords*sizeof(int)); recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
int *recLevelID = (int*) Malloc(nrecords*sizeof(int));
dtlist_type *dtlist = dtlist_new(); dtlist_type *dtlist = dtlist_new();
dtlist_set_stat(dtlist, timestat_date); dtlist_set_stat(dtlist, timestat_date);
...@@ -108,7 +113,7 @@ void *Timselstat(void *argument) ...@@ -108,7 +113,7 @@ void *Timselstat(void *argument)
field_type **samp1 = field_malloc(vlistID1, FIELD_NONE); field_type **samp1 = field_malloc(vlistID1, FIELD_NONE);
field_type **vars1 = field_malloc(vlistID1, FIELD_PTR); field_type **vars1 = field_malloc(vlistID1, FIELD_PTR);
field_type **vars2 = NULL; field_type **vars2 = NULL;
if ( lvarstd ) vars2 = field_malloc(vlistID1, FIELD_PTR); if ( lvarstd || lrange ) vars2 = field_malloc(vlistID1, FIELD_PTR);
for ( tsID = 0; tsID < noffset; tsID++ ) for ( tsID = 0; tsID < noffset; tsID++ )
{ {
...@@ -121,13 +126,13 @@ void *Timselstat(void *argument) ...@@ -121,13 +126,13 @@ void *Timselstat(void *argument)
if ( tsID == 0 ) if ( tsID == 0 )
{ {
recVarID[recID] = varID; recinfo[recID].varID = varID;
recLevelID[recID] = levelID; recinfo[recID].levelID = levelID;
} }
} }
} }
int otsID = 0; int otsID = 0;
if ( tsID < noffset ) if ( tsID < noffset )
{ {
cdoWarning("noffset is larger than number of timesteps!"); cdoWarning("noffset is larger than number of timesteps!");
...@@ -149,16 +154,25 @@ void *Timselstat(void *argument) ...@@ -149,16 +154,25 @@ void *Timselstat(void *argument)
if ( tsID == 0 ) if ( tsID == 0 )
{ {
recVarID[recID] = varID; recinfo[recID].varID = varID;
recLevelID[recID] = levelID; recinfo[recID].levelID = levelID;
} }
field_type *pvars1 = &vars1[varID][levelID];
field_type *pvars2 = vars2 ? &vars2[varID][levelID] : NULL;
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID)); gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
if ( nsets == 0 ) if ( nsets == 0 )
{ {
streamReadRecord(streamID1, vars1[varID][levelID].ptr, &nmiss); streamReadRecord(streamID1, pvars1->ptr, &nmiss);
vars1[varID][levelID].nmiss = (size_t)nmiss; pvars1->nmiss = (size_t)nmiss;
if ( lrange )
{
pvars2->nmiss = (size_t)nmiss;
for ( int i = 0; i < gridsize; i++ )
pvars2->ptr[i] = pvars1->ptr[i];
}
if ( nmiss > 0 || samp1[varID][levelID].ptr ) if ( nmiss > 0 || samp1[varID][levelID].ptr )
{ {
...@@ -166,19 +180,15 @@ void *Timselstat(void *argument) ...@@ -166,19 +180,15 @@ void *Timselstat(void *argument)
samp1[varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double)); samp1[varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
for ( int i = 0; i < gridsize; i++ ) for ( int i = 0; i < gridsize; i++ )
if ( DBL_IS_EQUAL(vars1[varID][levelID].ptr[i], samp1[varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
vars1[varID][levelID].missval) )
samp1[varID][levelID].ptr[i] = 0;
else
samp1[varID][levelID].ptr[i] = 1;
} }
} }
else else
{ {
streamReadRecord(streamID1, field.ptr, &nmiss); streamReadRecord(streamID1, field.ptr, &nmiss);
field.nmiss = (size_t)nmiss; field.nmiss = (size_t)nmiss;
field.grid = vars1[varID][levelID].grid; field.grid = pvars1->grid;
field.missval = vars1[varID][levelID].missval; field.missval = pvars1->missval;
if ( field.nmiss > 0 || samp1[varID][levelID].ptr ) if ( field.nmiss > 0 || samp1[varID][levelID].ptr )
{ {
...@@ -190,29 +200,38 @@ void *Timselstat(void *argument) ...@@ -190,29 +200,38 @@ void *Timselstat(void *argument)
} }
for ( int i = 0; i < gridsize; i++ ) for ( int i = 0; i < gridsize; i++ )
if ( !DBL_IS_EQUAL(field.ptr[i], vars1[varID][levelID].missval) ) if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
samp1[varID][levelID].ptr[i]++; samp1[varID][levelID].ptr[i]++;
} }
if ( lvarstd ) if ( lvarstd )
{ {
farsumq(&vars2[varID][levelID], field); farsumq(pvars2, field);
farsum(&vars1[varID][levelID], field); farsum(pvars1, field);
} }
else if ( lrange )
{
farmin(pvars2, field);
farmax(pvars1, field);
}
else else
{ {
farfun(&vars1[varID][levelID], field, operfunc); farfun(pvars1, field, operfunc);
} }
} }
} }
if ( nsets == 0 && lvarstd ) if ( nsets == 0 && lvarstd )
for ( varID = 0; varID < nvars; varID++ ) for ( int recID = 0; recID < maxrecs; recID++ )
{ {
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvars1 = &vars1[varID][levelID];
field_type *pvars2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT )