Skip to content

GitLab

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

New operator: merstd1 - Meridional standard deviation [Divisor is (n-1)]

parent fdbe9bf3
Hide whitespace changes
Inline Side-by-side
ChangeLog
View file @ 4e55e4a7
......@@ -5,6 +5,8 @@
2015-05-26 Uwe Schulzweida
* New operator: merstd1 - Meridional standard deviation [Divisor is (n-1)]
* New operator: mervar1 - Meridional variance [Divisor is (n-1)]
* New operator: zonstd1 - Zonal standard deviation [Divisor is (n-1)]
* New operator: zonvar1 - Zonal variance [Divisor is (n-1)]
......
OPERATORS
View file @ 4e55e4a7
......@@ -265,8 +265,10 @@ Operator catalog:
Merstat mersum Meridional sum
Merstat mermean Meridional mean
Merstat meravg Meridional average
Merstat mervar Meridional variance
Merstat merstd Meridional standard deviation
Merstat merstd1 Meridional standard deviation
Merstat mervar Meridional variance
Merstat mervar1 Meridional variance
Merstat merpctl Meridional percentiles
Gridboxstat gridboxmin Gridbox minimum
Gridboxstat gridboxmax Gridbox maximum
......
doc/tex/mod/Merstat
View file @ 4e55e4a7
......@@ -5,7 +5,7 @@
@Section = Statistical values
@Class = Statistic
@Arguments = ifile ofile
@Operators = mermin mermax mersum mermean meravg mervar merstd merpctl
@Operators = mermin mermax mersum mermean meravg merstd merstd1 mervar mervar1 merpctl
@BeginDescription
This module computes meridional statistical values of the input fields.
......@@ -65,7 +65,16 @@ For every longitude the area weighted average over all latitudes is computed.
@Title = Meridional variance
@BeginDescription
For every longitude the variance over all latitudes is computed.
For every longitude the variance over all latitudes is computed. Divisor is n.
@EndDescription
@EndOperator
@BeginOperator_mervar1
@Title = Meridional variance
@BeginDescription
For every longitude the variance over all latitudes is computed. Divisor is (n-1).
@EndDescription
@EndOperator
......@@ -74,7 +83,16 @@ For every longitude the variance over all latitudes is computed.
@Title = Meridional standard deviation
@BeginDescription
For every longitude the standard deviation over all latitudes is computed.
For every longitude the standard deviation over all latitudes is computed. Divisor is n.
@EndDescription
@EndOperator
@BeginOperator_merstd1
@Title = Meridional standard deviation
@BeginDescription
For every longitude the standard deviation over all latitudes is computed. Divisor is (n-1).
@EndDescription
@EndOperator
......
doc/tex/mod/Zonstat
View file @ 4e55e4a7
......@@ -83,7 +83,7 @@ For every latitude the variance over all longitudes is computed. Divisor is (n-1
@Title = Zonal standard deviation
@BeginDescription
For every latitude the standard deviation over all longitudes is computed. Divisor is n.
For every latitude the standard deviation over all longitudes is computed. Divisor is n.
@EndDescription
@EndOperator
......
src/Merstat.c
View file @ 4e55e4a7
......@@ -24,7 +24,9 @@
Merstat mermean Meridional mean
Merstat meravg Meridional average
Merstat merstd Meridional standard deviation
Merstat merstd Meridional standard deviation [Divisor is (n-1)]
Merstat mervar Meridional variance
Merstat mervar Meridional variance [Divisor is (n-1)]
Merstat merpctl Meridional percentiles
*/
......@@ -38,22 +40,13 @@
void *Merstat(void *argument)
{
int operatorID;
int operfunc;
int streamID1, streamID2;
int vlistID1, vlistID2;
int gridID1, gridID2 = -1, lastgrid = -1;
int wstatus = FALSE;
int nlonmax;
int index, ngrids;
int index;
int recID, nrecs;
int tsID, varID, levelID;
int lim;
int ndiffgrids;
int taxisID1, taxisID2;
int varID, levelID;
int needWeights = FALSE;
int pn = 0;
field_t field1, field2;
char varname[CDI_MAX_NAME];
cdoInitialize(argument);
......@@ -64,11 +57,13 @@ void *Merstat(void *argument)
cdoOperatorAdd("mermean", func_mean, 0, NULL);
cdoOperatorAdd("meravg", func_avg, 0, NULL);
cdoOperatorAdd("mervar", func_var, 0, NULL);
cdoOperatorAdd("mervar1", func_var1, 0, NULL);
cdoOperatorAdd("merstd", func_std, 0, NULL);
cdoOperatorAdd("merstd1", func_std1, 0, NULL);
cdoOperatorAdd("merpctl", func_pctl, 0, NULL);
operatorID = cdoOperatorID();
operfunc = cdoOperatorF1(operatorID);
int operatorID = cdoOperatorID();
int operfunc = cdoOperatorF1(operatorID);
/* RQ */
if ( operfunc == func_pctl )
......@@ -82,20 +77,21 @@ void *Merstat(void *argument)
/* QR */
if ( operfunc == func_mean || operfunc == func_avg ||
operfunc == func_var || operfunc == func_std )
operfunc == func_var || operfunc == func_std ||
operfunc == func_var1 || operfunc == func_std1 )
needWeights = TRUE;
streamID1 = streamOpenRead(cdoStreamName(0));
int streamID1 = streamOpenRead(cdoStreamName(0));
vlistID1 = streamInqVlist(streamID1);
vlistID2 = vlistDuplicate(vlistID1);
int vlistID1 = streamInqVlist(streamID1);
int vlistID2 = vlistDuplicate(vlistID1);
taxisID1 = vlistInqTaxis(vlistID1);
taxisID2 = taxisDuplicate(taxisID1);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
ngrids = vlistNgrids(vlistID1);
ndiffgrids = 0;
int ngrids = vlistNgrids(vlistID1);
int ndiffgrids = 0;
for ( index = 1; index < ngrids; index++ )
if ( vlistGrid(vlistID1, 0) != vlistGrid(vlistID1, index))
ndiffgrids++;
......@@ -118,17 +114,18 @@ void *Merstat(void *argument)
vlistChangeGridIndex(vlistID2, index, gridID2);
streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
streamDefVlist(streamID2, vlistID2);
gridID1 = vlistInqVarGrid(vlistID1, 0);
nlonmax = gridInqXsize(gridID1); /* max nlon ? */
int nlonmax = gridInqXsize(gridID1); /* max nlon ? */
int lim = vlistGridsizeMax(vlistID1);
field_t field1, field2;
field_init(&field1);
field_init(&field2);
lim = vlistGridsizeMax(vlistID1);
field1.ptr = (double*) malloc(lim*sizeof(double));
field1.weight = NULL;
if ( needWeights )
......@@ -137,7 +134,7 @@ void *Merstat(void *argument)
field2.ptr = (double*) malloc(nlonmax*sizeof(double));
field2.grid = gridID2;
tsID = 0;
int tsID = 0;
while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
{
taxisCopyTimestep(taxisID2, taxisID1);
......
src/field.c
View file @ 4e55e4a7
......@@ -332,7 +332,7 @@ double fldavg(field_t field)
static
void prevarsum(const double *restrict array, const double *restrict w, size_t len, int nmiss,
double missval, double *rsum, double *rsumw, double *rsumq, double *rsumwq)
double missval, double *restrict rsum, double *restrict rsumw, double *restrict rsumq, double *restrict rsumwq)
{
size_t i;
*rsum = *rsumw = 0;
......
src/field.h
View file @ 4e55e4a7
......@@ -117,7 +117,9 @@ void mersum(field_t field1, field_t *field2);
void meravg(field_t field1, field_t *field2);
void mermean(field_t field1, field_t *field2);
void merstd(field_t field1, field_t *field2);
void merstd1(field_t field1, field_t *field2);
void mervar(field_t field1, field_t *field2);
void mervar1(field_t field1, field_t *field2);
void merpctl(field_t field1, field_t *field2, const int k);
void fldrms(field_t field1, field_t field2, field_t *field3);
......
src/fieldmer.c
View file @ 4e55e4a7
......@@ -18,9 +18,8 @@
#include "cdo.h"
#include "cdo_int.h"
#include <cdi.h>
/* RQ */
#include "nth_element.h"
/* QR */
void merfun(field_t field1, field_t *field2, int function)
{
......@@ -30,7 +29,9 @@ void merfun(field_t field1, field_t *field2, int function)
else if ( function == func_mean ) mermean(field1, field2);
else if ( function == func_avg ) meravg(field1, field2);
else if ( function == func_std ) merstd(field1, field2);
else if ( function == func_std1 ) merstd1(field1, field2);
else if ( function == func_var ) mervar(field1, field2);
else if ( function == func_var1 ) mervar1(field1, field2);
else cdoAbort("function %d not implemented!", function);
}
......@@ -262,10 +263,42 @@ void meravg(field_t field1, field_t *field2)
field2->nmiss = rnmiss;
}
static
void prevarsum_mer(const double *restrict array, const double *restrict w, int nx, int ny, int nmiss,
double missval, double *restrict rsum, double *restrict rsumw, double *restrict rsumq, double *restrict rsumwq)
{
*rsum = 0;
*rsumq = 0;
*rsumw = 0;
*rsumwq = 0;
if ( nmiss > 0 )
{
for ( int j = 0; j < ny; j++ )
if ( !DBL_IS_EQUAL(array[j*nx], missval) &&
!DBL_IS_EQUAL(w[j*nx], missval) )
{
*rsum += w[j*nx] * array[j*nx];
*rsumq += w[j*nx] * array[j*nx] * array[j*nx];
*rsumw += w[j*nx];
*rsumwq += w[j*nx] * w[j*nx];
}
}
else
{
for ( int j = 0; j < ny; j++ )
{
*rsum += w[j*nx] * array[j*nx];
*rsumq += w[j*nx] * array[j*nx] * array[j*nx];
*rsumw += w[j*nx];
*rsumwq += w[j*nx] * w[j*nx];
}
}
}
void mervar(field_t field1, field_t *field2)
{
long i, j, nx, ny;
int rnmiss = 0;
int grid = field1.grid;
int nmiss = field1.nmiss;
......@@ -275,37 +308,12 @@ void mervar(field_t field1, field_t *field2)
double rsum = 0, rsumw = 0, rvar = 0;
double rsumq = 0, rsumwq = 0;
nx = gridInqXsize(grid);
ny = gridInqYsize(grid);
int nx = gridInqXsize(grid);
int ny = gridInqYsize(grid);
for ( i = 0; i < nx; i++ )
for ( int i = 0; i < nx; i++ )
{
rsum = 0;
rsumq = 0;
rsumw = 0;
rsumwq = 0;
if ( nmiss > 0 )
{
for ( j = 0; j < ny; j++ )
if ( !DBL_IS_EQUAL(array[j*nx+i], missval) &&
!DBL_IS_EQUAL(w[j*nx+i], missval) )
{
rsum += w[j*nx+i] * array[j*nx+i];
rsumq += w[j*nx+i] * array[j*nx+i] * array[j*nx+i];
rsumw += w[j*nx+i];
rsumwq += w[j*nx+i] * w[j*nx+i];
}
}
else
{
for ( j = 0; j < ny; j++ )
{
rsum += w[j*nx+i] * array[j*nx+i];
rsumq += w[j*nx+i] * array[j*nx+i] * array[j*nx+i];
rsumw += w[j*nx+i];
rsumwq += w[j*nx+i] * w[j*nx+i];
}
}
prevarsum_mer(array+i, w+i, nx, ny, nmiss, missval, &rsum, &rsumw, &rsumq, &rsumwq);
rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq*rsumw - rsum*rsum) / (rsumw*rsumw) : missval;
if ( rvar < 0 && rvar > -1.e-5 ) rvar = 0;
......@@ -319,32 +327,76 @@ void mervar(field_t field1, field_t *field2)
}
void mervar1(field_t field1, field_t *field2)
{
int rnmiss = 0;
int grid = field1.grid;
int nmiss = field1.nmiss;
double missval = field1.missval;
double *array = field1.ptr;
double *w = field1.weight;
double rsum = 0, rsumw = 0, rvar = 0;
double rsumq = 0, rsumwq = 0;
int nx = gridInqXsize(grid);
int ny = gridInqYsize(grid);
for ( int i = 0; i < nx; i++ )
{
prevarsum_mer(array+i, w+i, nx, ny, nmiss, missval, &rsum, &rsumw, &rsumq, &rsumwq);
rvar = (rsumw*rsumw > rsumwq) ? (rsumq*rsumw - rsum*rsum) / (rsumw*rsumw - rsumwq) : missval;
if ( rvar < 0 && rvar > -1.e-5 ) rvar = 0;
if ( DBL_IS_EQUAL(rvar, missval) ) rnmiss++;
field2->ptr[i] = rvar;
}
field2->nmiss = rnmiss;
}
void merstd(field_t field1, field_t *field2)
{
long i, nx;
int rnmiss = 0;
int grid = field1.grid;
double missval = field1.missval;
double rvar, rstd;
double rstd;
nx = gridInqXsize(grid);
int nx = gridInqXsize(grid);
mervar(field1, field2);
for ( i = 0; i < nx; i++ )
for ( int i = 0; i < nx; i++ )
{
rvar = field2->ptr[i];
rstd = var_to_std(field2->ptr[i], missval);
if ( DBL_IS_EQUAL(rvar, missval) || rvar < 0 )
{
rstd = missval;
}
else
{
rstd = IS_NOT_EQUAL(rvar, 0) ? sqrt(rvar) : 0;
}
if ( DBL_IS_EQUAL(rstd, missval) ) rnmiss++;
if ( DBL_IS_EQUAL(rvar, missval) ) rnmiss++;
field2->ptr[i] = rstd;
}
field2->nmiss = rnmiss;
}
void merstd1(field_t field1, field_t *field2)
{
int rnmiss = 0;
int grid = field1.grid;
double missval = field1.missval;
double rstd;
int nx = gridInqXsize(grid);
mervar1(field1, field2);
for ( int i = 0; i < nx; i++ )
{
rstd = var_to_std(field2->ptr[i], missval);
if ( DBL_IS_EQUAL(rstd, missval) ) rnmiss++;
field2->ptr[i] = rstd;
}
......
src/fieldzon.c
View file @ 4e55e4a7
......@@ -18,9 +18,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include <cdi.h>
/* RQ */
#include "nth_element.h"
/* QR */
void zonfun(field_t field1, field_t *field2, int function)
......@@ -365,10 +363,10 @@ void zonvar(field_t field1, field_t *field2)
double rsum = 0, rsumw = 0, rvar = 0;
double rsumq = 0, rsumwq = 0;
long nx = gridInqXsize(grid);
long ny = gridInqYsize(grid);
int nx = gridInqXsize(grid);
int ny = gridInqYsize(grid);
for ( long j = 0; j < ny; j++ )
for ( int j = 0; j < ny; j++ )
{
prevarsum_zon(array+j*nx, nx, nmiss, missval1, &rsum, &rsumw, &rsumq, &rsumwq);
......@@ -394,10 +392,10 @@ void zonvar1(field_t field1, field_t *field2)
double rsum = 0, rsumw = 0, rvar = 0;
double rsumq = 0, rsumwq = 0;
long nx = gridInqXsize(grid);
long ny = gridInqYsize(grid);
int nx = gridInqXsize(grid);
int ny = gridInqYsize(grid);
for ( long j = 0; j < ny; j++ )
for ( int j = 0; j < ny; j++ )
{
prevarsum_zon(array+j*nx, nx, nmiss, missval1, &rsum, &rsumw, &rsumq, &rsumwq);
......
src/modules.c
View file @ 4e55e4a7
......@@ -358,7 +358,7 @@ void *Maggraph(void *argument);
#define MergeOperators {"merge"}
#define MergegridOperators {"mergegrid"}
#define MergetimeOperators {"mergetime"}
#define MerstatOperators {"mermin", "mermax", "mersum", "mermean", "meravg", "mervar", "merstd", "merpctl"}
#define MerstatOperators {"mermin", "mermax", "mersum", "mermean", "meravg", "merstd", "merstd1", "mervar", "mervar1", "merpctl"}
#define MonarithOperators {"monadd", "monsub", "monmul", "mondiv"}
#define MrotuvOperators {"mrotuv"}
#define MrotuvbOperators {"mrotuvb"}
......
src/operator_help.h
View file @ 4e55e4a7
......@@ -1895,8 +1895,8 @@ static char *ZonstatHelp[] = {
static char *MerstatHelp[] = {
"NAME",
" mermin, mermax, mersum, mermean, meravg, mervar, merstd, merpctl - ",
" Meridional statistical values",
" mermin, mermax, mersum, mermean, meravg, merstd, merstd1, mervar, mervar1, ",
" merpctl - Meridional statistical values",
"",
"SYNOPSIS",
" <operator> ifile ofile",
......@@ -1919,10 +1919,14 @@ static char *MerstatHelp[] = {
" For every longitude the area weighted mean over all latitudes is computed.",
" meravg Meridional average",
" For every longitude the area weighted average over all latitudes is computed.",
" mervar Meridional variance",
" For every longitude the variance over all latitudes is computed.",
" merstd Meridional standard deviation",
" For every longitude the standard deviation over all latitudes is computed.",
" For every longitude the standard deviation over all latitudes is computed. Divisor is n.",
" merstd1 Meridional standard deviation",
" For every longitude the standard deviation over all latitudes is computed. Divisor is (n-1).",
" mervar Meridional variance",
" For every longitude the variance over all latitudes is computed. Divisor is n.",
" mervar1 Meridional variance",
" For every longitude the variance over all latitudes is computed. Divisor is (n-1).",
" merpctl Meridional percentiles",
" For every longitude the pth percentile over all latitudes is computed.",
"",
......
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

Imprint | Privacy Policy