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

Added new operator seasrange: seasonal range.

parent 87e4962a
......@@ -3,6 +3,10 @@
* Using CDI library version 1.9.0
* Version 1.9.0 release
2017-05-30 Uwe Schulzweida
* New operator seasrange: seasonal range (seasmax-seasmin)
2017-05-26 Uwe Schulzweida
* eca_gsl: set default to northern hemisphere
......
......@@ -18,6 +18,7 @@
/*
This module contains the following operators:
Seasstat seasrange Seasonal range
Seasstat seasmin Seasonal minimum
Seasstat seasmax Seasonal maximum
Seasstat seassum Seasonal sum
......@@ -34,7 +35,12 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "util.h"
typedef struct {
short varID;
short levelID;
} recinfo_t;
void *Seasstat(void *argument)
......@@ -46,22 +52,22 @@ void *Seasstat(void *argument)
int varID, levelID;
int year, month, day, seas0 = 0;
int nmiss;
int nlevel;
int oldmon = 0;
int nseason = 0;
const char *seas_name[4];
cdoInitialize(argument);
cdoOperatorAdd("seasmin", func_min, 0, NULL);
cdoOperatorAdd("seasmax", func_max, 0, NULL);
cdoOperatorAdd("seassum", func_sum, 0, NULL);
cdoOperatorAdd("seasmean", func_mean, 0, NULL);
cdoOperatorAdd("seasavg", func_avg, 0, NULL);
cdoOperatorAdd("seasvar", func_var, 0, NULL);
cdoOperatorAdd("seasvar1", func_var1, 0, NULL);
cdoOperatorAdd("seasstd", func_std, 0, NULL);
cdoOperatorAdd("seasstd1", func_std1, 0, NULL);
cdoOperatorAdd("seasrange", func_range, 0, NULL);
cdoOperatorAdd("seasmin", func_min, 0, NULL);
cdoOperatorAdd("seasmax", func_max, 0, NULL);
cdoOperatorAdd("seassum", func_sum, 0, NULL);
cdoOperatorAdd("seasmean", func_mean, 0, NULL);
cdoOperatorAdd("seasavg", func_avg, 0, NULL);
cdoOperatorAdd("seasvar", func_var, 0, NULL);
cdoOperatorAdd("seasvar1", func_var1, 0, NULL);
cdoOperatorAdd("seasstd", func_std, 0, NULL);
cdoOperatorAdd("seasstd1", func_std1, 0, NULL);
int operatorID = cdoOperatorID();
int operfunc = cdoOperatorF1(operatorID);
......@@ -69,10 +75,11 @@ void *Seasstat(void *argument)
int season_start = get_season_start();
get_season_name(seas_name);
int lmean = operfunc == func_mean || operfunc == func_avg;
int lstd = operfunc == func_std || operfunc == func_std1;
int lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
bool lrange = operfunc == func_range;
bool lmean = operfunc == func_mean || operfunc == func_avg;
bool lstd = operfunc == func_std || operfunc == func_std1;
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));
......@@ -88,11 +95,9 @@ void *Seasstat(void *argument)
streamDefVlist(streamID2, vlistID2);
int nvars = vlistNvars(vlistID1);
int nrecords = vlistNrecs(vlistID1);
int maxrecs = vlistNrecs(vlistID1);
int *recVarID = (int*) Malloc(nrecords*sizeof(int));
int *recLevelID = (int*) Malloc(nrecords*sizeof(int));
recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
dtlist_type *dtlist = dtlist_new();
dtlist_set_stat(dtlist, timestat_date);
......@@ -107,10 +112,10 @@ void *Seasstat(void *argument)
field_type **samp1 = field_malloc(vlistID1, FIELD_NONE);
field_type **vars1 = field_malloc(vlistID1, FIELD_PTR);
field_type **vars2 = NULL;
if ( lvarstd ) vars2 = field_malloc(vlistID1, FIELD_PTR);
if ( lvarstd || lrange ) vars2 = field_malloc(vlistID1, FIELD_PTR);
int tsID = 0;
int otsID = 0;
int tsID = 0;
int otsID = 0;
while ( TRUE )
{
long nsets = 0;
......@@ -150,16 +155,25 @@ void *Seasstat(void *argument)
if ( tsID == 0 )
{
recVarID[recID] = varID;
recLevelID[recID] = levelID;
recinfo[recID].varID = varID;
recinfo[recID].levelID = levelID;
}
field_type *pvar1 = &vars1[varID][levelID];
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
if ( nsets == 0 )
{
streamReadRecord(streamID1, vars1[varID][levelID].ptr, &nmiss);
vars1[varID][levelID].nmiss = nmiss;
streamReadRecord(streamID1, pvar1->ptr, &nmiss);
pvar1->nmiss = (size_t)nmiss;
if ( lrange )
{
field_type *pvar2 = &vars2[varID][levelID];
for ( int i = 0; i < gridsize; i++ )
pvar2->ptr[i] = pvar1->ptr[i];
pvar2->nmiss = (size_t)nmiss;
}
if ( nmiss > 0 || samp1[varID][levelID].ptr )
{
......@@ -167,19 +181,15 @@ void *Seasstat(void *argument)
samp1[varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
for ( int i = 0; i < gridsize; i++ )
if ( DBL_IS_EQUAL(vars1[varID][levelID].ptr[i],
vars1[varID][levelID].missval) )
samp1[varID][levelID].ptr[i] = 0;
else
samp1[varID][levelID].ptr[i] = 1;
samp1[varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvar1->ptr[i], pvar1->missval);
}
}
else
{
streamReadRecord(streamID1, field.ptr, &nmiss);
field.nmiss = (size_t)nmiss;
field.grid = vars1[varID][levelID].grid;
field.missval = vars1[varID][levelID].missval;
field.grid = pvar1->grid;
field.missval = pvar1->missval;
if ( field.nmiss > 0 || samp1[varID][levelID].ptr )
{
......@@ -191,29 +201,40 @@ void *Seasstat(void *argument)
}
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], pvar1->missval) )
samp1[varID][levelID].ptr[i]++;
}
if ( lvarstd )
{
farsumq(&vars2[varID][levelID], field);
farsum(&vars1[varID][levelID], field);
field_type *pvar2 = &vars2[varID][levelID];
farsumq(pvar2, field);
farsum(pvar1, field);
}
else if ( lrange )
{
field_type *pvar2 = &vars2[varID][levelID];
farmin(pvar2, field);
farmax(pvar1, field);
}
else
{
farfun(&vars1[varID][levelID], field, operfunc);
farfun(pvar1, field, operfunc);
}
}
}
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 *pvar1 = &vars1[varID][levelID];
field_type *pvar2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
farmoq(&vars2[varID][levelID], vars1[varID][levelID]);
farmoq(pvar2, *pvar1);
}
vdate1 = vdate;
......@@ -225,41 +246,52 @@ void *Seasstat(void *argument)
if ( nrecs == 0 && nsets == 0 ) break;
if ( lmean )
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 *pvar1 = &vars1[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
if ( samp1[varID][levelID].ptr == NULL )
farcdiv(&vars1[varID][levelID], (double)nsets);
else
fardiv(&vars1[varID][levelID], samp1[varID][levelID]);
}
if ( samp1[varID][levelID].ptr == NULL )
farcdiv(pvar1, (double)nsets);
else
fardiv(pvar1, samp1[varID][levelID]);
}
else if ( lvarstd )
for ( varID = 0; varID < nvars; varID++ )
{
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
if ( samp1[varID][levelID].ptr == NULL )
{
if ( lstd )
farcstd(&vars1[varID][levelID], vars2[varID][levelID], nsets, divisor);
else
farcvar(&vars1[varID][levelID], vars2[varID][levelID], nsets, divisor);
}
else
{
if ( lstd )
farstd(&vars1[varID][levelID], vars2[varID][levelID], samp1[varID][levelID], divisor);
else
farvar(&vars1[varID][levelID], vars2[varID][levelID], samp1[varID][levelID], divisor);
}
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvar1 = &vars1[varID][levelID];
field_type *pvar2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
if ( samp1[varID][levelID].ptr == NULL )
{
if ( lstd ) farcstd(pvar1, *pvar2, nsets, divisor);
else farcvar(pvar1, *pvar2, nsets, divisor);
}
else
{
if ( lstd ) farstd(pvar1, *pvar2, samp1[varID][levelID], divisor);
else farvar(pvar1, *pvar2, samp1[varID][levelID], divisor);
}
}
else if ( lrange )
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvar1 = &vars1[varID][levelID];
field_type *pvar2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
farsub(pvar1, *pvar2);
}
if ( cdoVerbose )
{
......@@ -284,15 +316,16 @@ void *Seasstat(void *argument)
otsID+1, vdatestr, nsets, nsets == 1 ? "" : "s");
}
for ( int recID = 0; recID < nrecords; recID++ )
for ( int recID = 0; recID < maxrecs; recID++ )
{
varID = recVarID[recID];
levelID = recLevelID[recID];
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvar1 = &vars1[varID][levelID];
if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, vars1[varID][levelID].ptr, (int)vars1[varID][levelID].nmiss);
streamWriteRecord(streamID2, pvar1->ptr, (int)pvar1->nmiss);
}
if ( nrecs == 0 ) break;
......@@ -306,8 +339,7 @@ void *Seasstat(void *argument)
if ( field.ptr ) Free(field.ptr);
if ( recVarID ) Free(recVarID);
if ( recLevelID ) Free(recLevelID);
Free(recinfo);
dtlist_delete(dtlist);
......
......@@ -160,7 +160,7 @@ void *Timstat(void *argument)
bool lmean = operfunc == func_mean || operfunc == func_avg;
bool lstd = operfunc == func_std || operfunc == func_std1;
bool lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
if ( operfunc == func_mean )
{
......@@ -311,10 +311,7 @@ void *Timstat(void *argument)
samp1[varID][levelID].ptr = (double*) Malloc(nwpv*gridsize*sizeof(double));
for ( int i = 0; i < nwpv*gridsize; i++ )
if ( DBL_IS_EQUAL(pvar1->ptr[i], pvar1->missval) )
samp1[varID][levelID].ptr[i] = 0;
else
samp1[varID][levelID].ptr[i] = 1;
samp1[varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvar1->ptr[i], pvar1->missval);
}
}
else
......
......@@ -406,7 +406,7 @@ void *Samplegrid(void *argument); // "samplegrid", "subgrid"
#define SamplegridiconOperators {"samplegridicon"}
#define SeascountOperators {"seascount"}
#define SeaspctlOperators {"seaspctl"}
#define SeasstatOperators {"seasmin", "seasmax", "seassum", "seasmean", "seasavg", "seasstd", "seasstd1", "seasvar", "seasvar1"}
#define SeasstatOperators {"seasrange", "seasmin", "seasmax", "seassum", "seasmean", "seasavg", "seasstd", "seasstd1", "seasvar", "seasvar1"}
#define SelboxOperators {"sellonlatbox", "selindexbox"}
#define SelgridcellOperators {"selgridcell", "delgridcell"}
#define SelectOperators {"select", "delete"}
......@@ -510,7 +510,7 @@ void *Samplegrid(void *argument); // "samplegrid", "subgrid"
#define YmonstatOperators {"ymonmin", "ymonmax", "ymonsum", "ymonmean", "ymonavg", "ymonstd", "ymonstd1", "ymonvar", "ymonvar1"}
#define YseaspctlOperators {"yseaspctl"}
#define YseasstatOperators {"yseasmin", "yseasmax", "yseassum", "yseasmean", "yseasavg", "yseasstd", "yseasstd1", "yseasvar", "yseasvar1"}
#define ZonstatOperators {"zonmin", "zonmax", "zonrange", "zonsum", "zonmean", "zonavg", "zonstd", "zonstd1", "zonvar", "zonvar1", "zonpctl"}
#define ZonstatOperators {"zonrange", "zonmin", "zonmax", "zonsum", "zonmean", "zonavg", "zonstd", "zonstd1", "zonvar", "zonvar1", "zonpctl"}
#define EcaCfdOperators {"eca_cfd"}
#define EcaCsuOperators {"eca_csu"}
......
#! @BASH@
echo 1..9 # Number of tests to be executed.
echo 1..10 # Number of tests to be executed.
#
test -n "$CDO" || CDO=cdo
test -n "$DATAPATH" || DATAPATH=./data
......@@ -7,7 +7,7 @@ test -n "$DATAPATH" || DATAPATH=./data
CDOOUT=cout$$
CDOERR=cerr$$
TYPE=seas
STATS="min max sum avg mean std std1 var var1"
STATS="min max range sum avg mean std std1 var var1"
NTEST=1
#
IFILE=$DATAPATH/ts_mm_5years
......
......@@ -9,7 +9,7 @@ COMP_REF = comp_eqc_ref comp_gec_ref comp_gtc_ref comp_lec_ref comp_ltc_ref
EOF_REF = eval_ref eof_ref pcoeff00000
YMONSTAT_REF = ymonmin_ref ymonmax_ref ymonsum_ref ymonavg_ref ymonmean_ref ymonstd_ref ymonstd1_ref ymonvar_ref ymonvar1_ref
YSEASSTAT_REF = yseasmin_ref yseasmax_ref yseassum_ref yseasavg_ref yseasmean_ref yseasstd_ref yseasstd1_ref yseasvar_ref yseasvar1_ref
SEASSTAT_REF = seasmin_ref seasmax_ref seassum_ref seasavg_ref seasmean_ref seasstd_ref seasstd1_ref seasvar_ref seasvar1_ref
SEASSTAT_REF = seasmin_ref seasmax_ref seassum_ref seasavg_ref seasmean_ref seasstd_ref seasstd1_ref seasvar_ref seasvar1_ref seasrange_ref
RUNSTAT_REF = runmin_ref runmax_ref runsum_ref runavg_ref runmean_ref runstd_ref runstd1_ref runvar_ref runvar1_ref
TIMSTAT_REF = timmin_ref timmax_ref timsum_ref timavg_ref timmean_ref timstd_ref timstd1_ref timvar_ref timvar1_ref timrange_ref
YEARSTAT_REF = yearmin_ref yearmax_ref yearsum_ref yearavg_ref yearmean_ref yearstd_ref yearstd1_ref yearvar_ref yearvar1_ref yearrange_ref
......
......@@ -299,7 +299,7 @@ COMP_REF = comp_eqc_ref comp_gec_ref comp_gtc_ref comp_lec_ref comp_ltc_ref comp
EOF_REF = eval_ref eof_ref pcoeff00000
YMONSTAT_REF = ymonmin_ref ymonmax_ref ymonsum_ref ymonavg_ref ymonmean_ref ymonstd_ref ymonstd1_ref ymonvar_ref ymonvar1_ref
YSEASSTAT_REF = yseasmin_ref yseasmax_ref yseassum_ref yseasavg_ref yseasmean_ref yseasstd_ref yseasstd1_ref yseasvar_ref yseasvar1_ref
SEASSTAT_REF = seasmin_ref seasmax_ref seassum_ref seasavg_ref seasmean_ref seasstd_ref seasstd1_ref seasvar_ref seasvar1_ref
SEASSTAT_REF = seasmin_ref seasmax_ref seassum_ref seasavg_ref seasmean_ref seasstd_ref seasstd1_ref seasvar_ref seasvar1_ref seasrange_ref
RUNSTAT_REF = runmin_ref runmax_ref runsum_ref runavg_ref runmean_ref runstd_ref runstd1_ref runvar_ref runvar1_ref
TIMSTAT_REF = timmin_ref timmax_ref timsum_ref timavg_ref timmean_ref timstd_ref timstd1_ref timvar_ref timvar1_ref timrange_ref
YEARSTAT_REF = yearmin_ref yearmax_ref yearsum_ref yearavg_ref yearmean_ref yearstd_ref yearstd1_ref yearvar_ref yearvar1_ref yearrange_ref
......
......@@ -25,9 +25,9 @@ STATS="min max sum avg mean std std1 var var1 range"
STATS="range"
#
IFILE=ts_mm_5years
#for STAT in $STATS; do
# $CDO $FORMAT seas${STAT} $IFILE seas${STAT}_ref
#done
for STAT in $STATS; do
$CDO $FORMAT seas${STAT} $IFILE seas${STAT}_ref
done
#
IFILE=ts_mm_5years
#for STAT in $STATS; do
......
Markdown is supported
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