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

Added operator eca_pd - Precipitation days index per time period

parent 98e85119
......@@ -90,6 +90,7 @@ doc/tex/mod/EcaHd -text
doc/tex/mod/EcaHwdi -text
doc/tex/mod/EcaHwfi -text
doc/tex/mod/EcaId -text
doc/tex/mod/EcaPd -text
doc/tex/mod/EcaR10mm -text
doc/tex/mod/EcaR20mm -text
doc/tex/mod/EcaR75p -text
......
2010-11-16 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* Added operator eca_pd - Precipitation days index per time period [request: Barbara Hennemuth]
2010-11-13 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* sellonlatbox: fix rounding error of the last lon index [report: Mondher Chekki]
......
......@@ -7,9 +7,10 @@ Version 1.4.7 (?? December 2010):
* improved support for netCDF output from WRF model (import time and grid variables)
New operators:
* ydayadd, ydaysub, ydaymul, ydaydiv: Multi-year daily arithmetic
* eca_pd: Precipitation days index per time period
* dv2ps: Divergence and vorticity to velocity potential and stream function
Fixed bugs:
* fldmean: gives wrong result for grid units radian
* fldmean: gives wrong result for grid units [radian]
* Yseasstat: bug fix for datasets with time constant fields
* sellevel: fix problem with hybrid model levels and netCDF output
* sellonlatbox: fix rounding error of the last lon index
......
......@@ -206,6 +206,8 @@ Operator catalog:
Ensstat enspctl Ensemble percentiles
Ensstat2 ensbrs Brier score
Ensstat2 enscrps Cumulative Ranked Probability score
Ensstat2 ensrkhistspace Ranked Histogram averaged over time
Ensstat2 ensrkhisttime Ranked Histogram averaged over space
Fldstat fldmin Field minimum
Fldstat fldmax Field maximum
Fldstat fldsum Field sum
......@@ -465,8 +467,9 @@ Operator catalog:
EcaHwdi eca_hwdi Heat wave duration index wrt mean of reference period
EcaHwfi eca_hwfi Warm spell days index wrt 90th percentile of reference period
EcaId eca_id Ice days index per time period
EcaR10mm eca_r10mm Heavy precipitation days index per time period
EcaR20mm eca_r20mm Very heavy precipitation days index per time period
EcaPd eca_pd Precipitation days index per time period
EcaPd eca_r10mm Heavy precipitation days index per time period
EcaPd eca_r20mm Very heavy precipitation days index per time period
EcaR75p eca_r75p Moderate wet days wrt 75th percentile of reference period
EcaR75ptot eca_r75ptot Precipitation percent due to R75p days
EcaR90p eca_r90p Wet days wrt 90th percentile of reference period
......
......@@ -125,8 +125,7 @@ EcaHd Climate indices
EcaHwdi Climate indices
EcaHwfi Climate indices
EcaId Climate indices
EcaR10mm Climate indices
EcaR20mm Climate indices
EcaPd Climate indices
EcaR75p Climate indices
EcaR75ptot Climate indices
EcaR90p Climate indices
......
......@@ -373,6 +373,7 @@ while (<MOFILE>) {
if ( $len > $maxlen ) {
$maxlen = $len;
$maxitem = $operator;
$maxitem =~ s/_/\\_/g;
}
}
$maxlen += 2;
......@@ -611,14 +612,14 @@ while (<MOFILE>) {
$otitle = @oper_title[$opercnt-1];
$oparameter = @oper_parameter[$opercnt-1];
$opernamex = $opername;
$opernamex =~ s/_/\\_/g;
if ( $#hkeys > 0 ) {
printf HELPFILE (" \" %-*s%s\",\n",
$maxlen, $operator, $otitle);
if ( $xopercnt == 1 ) {
print TRFILE "\\begin{defalist}{\\bf $maxitem \\ }\n";
}
$operatorx = $operator;
$operatorx =~ s/_/\\_/g;
print TRFILE "\\item[{\\bf $operatorx}\\ \\ \\hfill]\n";
print TRFILE "$otitle \\\\\n";
}
......@@ -628,7 +629,7 @@ while (<MOFILE>) {
print TRCARD "\\rowcolor[gray]{.9}\n";
print TRCARD "\\makebox[$len2][r]{Syntax} ";
if ( $xopercnt == 2 ) {
print TRCARD "& \\hspace*{0mm}{\\bf $opername}{\\sl $operpara} \ {\\tt $marguments} ";
print TRCARD "& \\hspace*{0mm}{\\bf $opernamex}{\\sl $operpara} \ {\\tt $marguments} ";
} else {
print TRCARD "& \\hspace*{0mm}{\$<\\!operator\\!>\$}{\\sl $operpara} \ {\\tt $marguments} ";
}
......@@ -751,6 +752,7 @@ while (<MOFILE>) {
if ( $len > $maxlen ) {
$maxlen = $len;
$maxitem = $value;
$maxitem =~ s/_/\\_/g;
}
}
}
......@@ -1002,7 +1004,7 @@ while (<MOFILE>) {
print TRCARD "& \\hspace*{0mm}{\$<\\!operator\\!>\$}{\\sl $operpara} \ {\\tt $marguments} ";
}
print TRCARD "\\\\ \\hline \n";
print TRCARD "\\end{tabular*}\n";
print TRCARD "\\vspace{1mm}\n";
print TRCARD "\n";
......
@BeginModule
@NewPage
@Name = EcaPd
@Title = Precipitation days index per time period
@Section = Climate indices
@Class = Climate index
@Arguments = ifile ofile
@Operators = eca_pd eca_r10mm eca_r20mm
@BeginDescription
Let @file{ifile} be a time series of daily precipitation amounts RR in [mm]
(or alternatively in [kg m-2]), then counted is the number of days where RR is at least @math{x} mm.
eca\_r10mm and eca\_r20mm are specific ECA operators with a daily precipitation amount of 10
and 20 mm respectively.
The date information of a time step in ofile is the date of the last
contributing time step in @file{ifile}.
@EndDescription
@EndModule
@BeginOperator_eca_pd
@Title = Precipitation days index per time period
@Parameter = x
@BeginDescription
Generic ECA operator with a daily precipitation sum exceeding @math{x} mm.
@EndDescription
@EndOperator
@BeginOperator_eca_r10mm
@Title = Heavy precipitation days index per time period
@BeginDescription
Specific ECA operator with a daily precipitation sum exceeding 10 mm.
@EndDescription
@EndOperator
@BeginOperator_eca_r20mm
@Title = Very heavy precipitation days index per time period
@BeginDescription
Specific ECA operator with a daily precipitation sum exceeding 20 mm.
@EndDescription
@EndOperator
@BeginParameter x
@Item = x
FLOAT Daily precipitation amount threshold in [mm]
@EndParameter
@BeginNote
Please note, that precipitation rates in [mm/s] or [kg m-2/s] have to be converted to
precipitation amounts (multiply with 86400 s).
@EndNote
@BeginExample
To get the number of days with precipitation greater than 25 mm for
a time series of daily precipitation amounts use:
@BeginVerbatim
cdo eca_pd,25 ifile ofile
@EndVerbatim
@EndExample
......@@ -176,9 +176,9 @@ static const char CWD_NAME2[] = "number_of_cwd_periods_with_more_than_5da
static const char CWD_LONGNAME2[] = "Number of cwd periods in given time period with more than 5 days. The time period should be defined by the bounds of the time coordinate.";
static const char CWD_UNITS2[] = "No.";
static const char RINDEX_NAME[] = "precipitation_days_index_per_time_period";
static const char RINDEX_LONGNAME[] = "precipitation days is the number of days per time period with daily precipitation sum exceeding RR mm. The time period should be defined by the bounds of the time coordinate.";
static const char RINDEX_UNITS[] = "No.";
static const char PD_NAME[] = "precipitation_days_index_per_time_period";
static const char PD_LONGNAME[] = "precipitation days is the number of days per time period with daily precipitation sum exceeding %g mm. The time period should be defined by the bounds of the time coordinate.";
static const char PD_UNITS[] = "No.";
static const char R10MM_NAME[] = "heavy_precipitation_days_index_per_time_period";
static const char R10MM_LONGNAME[] = "Heavy precipitation days is the number of days per time period with daily precipitation sum exceeding 10 mm. The time period should be defined by the bounds of the time coordinate.";
......@@ -942,39 +942,41 @@ void *EcaCwd(void *argument)
}
void *EcaRainIndex(void *argument)
void *EcaPd(void *argument)
{
int ECA_RAININDEX, ECA_R10MM, ECA_R20MM;
int ECA_PD, ECA_R10MM, ECA_R20MM;
int operatorID;
double rr = 0;
char lnamebuffer[1024];
double threshold = 0;
ECA_REQUEST_1 request;
cdoInitialize(argument);
ECA_RAININDEX = cdoOperatorAdd("eca_rainindex", 0, 31, NULL);
ECA_R10MM = cdoOperatorAdd("eca_r10mm", 0, 31, NULL);
ECA_R20MM = cdoOperatorAdd("eca_r20mm", 0, 31, NULL);
ECA_PD = cdoOperatorAdd("eca_pd", 0, 31, NULL);
ECA_R10MM = cdoOperatorAdd("eca_r10mm", 0, 31, NULL);
ECA_R20MM = cdoOperatorAdd("eca_r20mm", 0, 31, NULL);
operatorID = cdoOperatorID();
if ( operatorID == ECA_RAININDEX )
if ( operatorID == ECA_PD )
{
operatorInputArg("RR");
operatorInputArg("daily precipitation amount threshold in [mm]");
if ( operatorArgc() < 1 ) cdoAbort("Not enough arguments!");
if ( operatorArgc() > 1 ) cdoAbort("Too many arguments!");
rr = atof(operatorArgv()[0]);
threshold = atof(operatorArgv()[0]);
if ( rr < 0 ) cdoAbort("Parameter out of range: RR = %d", rr);
if ( threshold < 0 ) cdoAbort("Parameter out of range: threshold = %d", threshold);
request.var1.name = RINDEX_NAME;
request.var1.longname = RINDEX_LONGNAME;
request.var1.units = RINDEX_UNITS;
sprintf(lnamebuffer, PD_LONGNAME, threshold);
request.var1.name = PD_NAME;
request.var1.longname = lnamebuffer;
request.var1.units = PD_UNITS;
}
else if ( operatorID == ECA_R10MM )
{
rr = 10;
threshold = 10;
request.var1.name = R10MM_NAME;
request.var1.longname = R10MM_LONGNAME;
......@@ -982,17 +984,17 @@ void *EcaRainIndex(void *argument)
}
else if ( operatorID == ECA_R20MM )
{
rr = 20;
threshold = 20;
request.var1.name = R20MM_NAME;
request.var1.longname = R20MM_LONGNAME;
request.var1.units = R20MM_UNITS;
}
if ( cdoVerbose ) cdoPrint("RR = %g", rr);
if ( cdoVerbose ) cdoPrint("threshold = %g", threshold);
request.var1.f1 = farselgec;
request.var1.f1arg = rr;
request.var1.f1arg = threshold;
request.var1.f2 = farnum;
request.var1.f3 = NULL;
request.var1.mulc = 0.0;
......
......@@ -220,7 +220,7 @@ void *EcaTx90p(void *argument);
void *EcaCdd(void *argument);
void *EcaCwd(void *argument);
void *EcaRr1(void *argument);
void *EcaRainIndex(void *argument);
void *EcaPd(void *argument);
void *EcaR75p(void *argument);
void *EcaR75ptot(void *argument);
void *EcaR90p(void *argument);
......@@ -455,9 +455,7 @@ void *Wct(void *argument);
#define EcaCddOperators {"eca_cdd"}
#define EcaCwdOperators {"eca_cwd"}
#define EcaRr1Operators {"eca_rr1"}
#define EcaRainIndexOperators {"eca_rainindex"}
#define EcaR10mmOperators {"eca_r10mm"}
#define EcaR20mmOperators {"eca_r20mm"}
#define EcaPdOperators {"eca_pd", "eca_r10mm", "eca_r20mm"}
#define EcaR75pOperators {"eca_r75p"}
#define EcaR75ptotOperators {"eca_r75ptot"}
#define EcaR90pOperators {"eca_r90p"}
......@@ -685,9 +683,7 @@ static modules_t Modules[] =
{ EcaCdd, EcaCddHelp, EcaCddOperators, CDI_REAL, 1, 1 },
{ EcaCwd, EcaCwdHelp, EcaCwdOperators, CDI_REAL, 1, 1 },
{ EcaRr1, EcaRr1Help, EcaRr1Operators, CDI_REAL, 1, 1 },
{ EcaRainIndex, NULL, EcaRainIndexOperators, CDI_REAL, 1, 1 },
{ EcaRainIndex, EcaR10mmHelp, EcaR10mmOperators, CDI_REAL, 1, 1 },
{ EcaRainIndex, EcaR20mmHelp, EcaR20mmOperators, CDI_REAL, 1, 1 },
{ EcaPd, EcaPdHelp, EcaPdOperators, CDI_REAL, 1, 1 },
{ EcaR75p, EcaR75pHelp, EcaR75pOperators, CDI_REAL, 2, 1 },
{ EcaR75ptot, EcaR75ptotHelp, EcaR75ptotOperators, CDI_REAL, 2, 1 },
{ EcaR90p, EcaR90pHelp, EcaR90pOperators, CDI_REAL, 2, 1 },
......
......@@ -1353,17 +1353,20 @@ static char *EnsstatHelp[] = {
static char *Ensstat2Help[] = {
"NAME",
" ensbrs, enscrps - Statistical values over an ensemble",
" ensbrs, enscrps, ensrkhistspace, ensrkhisttime - ",
" Statistical values over an ensemble",
"",
"SYNOPSIS",
" <operator> obsfile ensfiles ofile",
"",
"DESCRIPTION",
" This module computes statistical values over the ensemble of ensfiles using",
" obsfile as a reference. Depending on the chosen operator the Brier score or ",
" Cumulative ranked probability score over all input files or a skill score is written ",
" to ofile. The date information of a time step in ofile is the date of the ",
" first input file.",
" obsfile as a reference. Depending on the chosen operator the Brier score, ",
" Cumulative ranked probability score or a ranked Histogram over all Ensembles ensfiles",
" with reference to obsfile is written to ofile. ",
" The date and grid information of a time step in ofile is the date of the ",
" first input file. Thus all input files are required to have the same structure in ",
" terms of the gridsize, variable definitions and number of time steps. ",
" ",
" All Operators in this module use obsfile as the reference (for instance ",
" an observation) whereas ensfiles are understood as an ensemble consisting ",
......@@ -1379,12 +1382,27 @@ static char *Ensstat2Help[] = {
" The cumulative ranked probability score is understood as the integral of the ",
" squared difference between the cumulative distribution function cdf(ensfiles(t,x)) ",
" of the ensemble and a heavyside function H(obsfile(t,x)). ",
" ",
" The operators ensrkhistspace and ensrkhisttime compute Ranked Histograms. ",
" Therefor the vertical axis is utilized as the Histogram axis, which prohibits",
" the use of files containing more than one level. The histogram axis has ",
" nensfiles+1 bins with level 0 containing for each grid point the number of ",
" observations being smaller as all ensembles and level nensfiles+1 indicating",
" the number of observations being larger than all ensembles. ",
" ",
" ensrkhistspace computes a ranked histogram at each time step reducing each ",
" horizontal grid to a 1x1 grid and keeping the time axis as in obsfile. ",
" Contrary ensrkhistspace computes a histogram at each grid point keeping the ",
" horizontal grid for each variable and reducing the time-axis. The time infor-",
" mation is that from the last time step in obsfile. ",
"",
"OPERATORS",
" ensbrs Brier score",
" o(t,x) = (i2(t,x)-i1(t,x))^2 + (i3(t,x)-i1(t,x))^2 + ... + (in(t,x)-i1(t,x))^2",
" enscrps Cumulative Ranked Probability score",
" o(t,x) = crps( {i1(t,x)}, {i2(t,x),...,in(t,x)} )",
" ensbrs Brier score",
" o(t,x) = (i2(t,x)-i1(t,x))^2 + (i3(t,x)-i1(t,x))^2 + ... + (in(t,x)-i1(t,x))^2",
" enscrps Cumulative Ranked Probability score",
" o(t,x) = crps( {i1(t,x)}, {i2(t,x),...,in(t,x)} )",
" ensrkhistspace Ranked Histogram averaged over time",
" ensrkhisttime Ranked Histogram averaged over space",
NULL
};
......@@ -4038,37 +4056,35 @@ static char *EcaIdHelp[] = {
NULL
};
static char *EcaR10mmHelp[] = {
static char *EcaPdHelp[] = {
"NAME",
" eca_r10mm - Heavy precipitation days index per time period",
" eca_pd, eca_r10mm, eca_r20mm - Precipitation days index per time period",
"",
"SYNOPSIS",
" eca_pd,x ifile ofile",
" eca_r10mm ifile ofile",
" eca_r20mm ifile ofile",
"",
"DESCRIPTION",
" Let ifile be a time series of daily precipitation amounts RR,",
" then counted is the number of days where RR is at least 10 mm. ",
" The date information of a time step in ofile is the date of",
" the last contributing time step in ifile.",
" Let ifile be a time series of daily precipitation amounts RR in [mm]",
" (or alternatively in [kg m-2]), then counted is the number of days where RR is at least x mm. ",
" eca\\_r10mm and eca\\_r20mm are specific ECA operators with a daily precipitation amount of 10 ",
" and 20 mm respectively.",
" The date information of a time step in ofile is the date of the last",
" contributing time step in ifile.",
" The following variables are created: ",
" - heavy_precipitation_days_index_per_time_period",
NULL
};
static char *EcaR20mmHelp[] = {
"NAME",
" eca_r20mm - Very heavy precipitation days index per time period",
" - precipitation_days_index_per_time_period",
"",
"SYNOPSIS",
" eca_r20mm ifile ofile",
"OPERATORS",
" eca_pd Precipitation days index per time period",
" Generic ECA operator with a daily precipitation sum exceeding x mm.",
" eca_r10mm Heavy precipitation days index per time period",
" Specific ECA operator with a daily precipitation sum exceeding 10 mm.",
" eca_r20mm Very heavy precipitation days index per time period",
" Specific ECA operator with a daily precipitation sum exceeding 20 mm.",
"",
"DESCRIPTION",
" Let ifile be a time series of daily precipitation amounts RR,",
" then counted is the number of days where RR is at least 20 mm. ",
" The date information of a time step in ofile is the date of",
" the last contributing time step in ifile.",
" The following variables are created: ",
" - very_heavy_precipitation_days_index_per_time_period",
"PARAMETER",
" x FLOAT Daily precipitation amount threshold in [mm]",
NULL
};
......
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