UFO
ProbabilityOfGrossError.cc
Go to the documentation of this file.
1 /*
2  * (C) Crown copyright 2020, Met Office
3  *
4  * This software is licensed under the terms of the Apache Licence Version 2.0
5  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6  */
7 
9 
10 namespace ufo {
12  const std::vector<float> &obsVal,
13  const std::vector<float> &obsErr,
14  const std::vector<float> &bkgVal,
15  const std::vector<float> &bkgErr,
16  const std::vector<float> &PdBad,
17  const bool ModelLevels,
18  std::vector<int> &flags,
19  std::vector<float> &PGE,
20  std::vector<float> &PGEBd,
21  float ErrVarMax,
22  const std::vector<float> *obsVal2,
23  const std::vector<float> *bkgVal2)
24  {
25  const float missingValueFloat = util::missingValue(1.0f);
26  // PGE multiplication factor used to store PGE values for later use.
27  const double PGEMult = 1000.0;
28  // Missing data indicator for stored PGEs.
29  const double PGEMDI = 1.111;
30  // Maximum value of exponent in background QC.
31  const double ExpArgMax = options.PGE_ExpArgMax.value();
32  // PGE rejection limit.
33  const double PGECrit = options.PGE_PGECrit.value();
34  // Multiplication factor for observation errors.
35  const float ObErrMult = options.PGE_ObErrMult.value();
36  // Multiplication factor for background errors.
37  const float BkgErrMult = options.PGE_BkgErrMult.value();
38  // Critical value for squared difference from background / ErrVar.
39  const double SDiffCrit = obsVal2 && bkgVal2 ?
40  options.PGE_SDiffCrit.value() * 2.0 :
41  options.PGE_SDiffCrit.value();
42  // Number of levels in profile, or total number of single-level obs.
43  const size_t numLocs = obsVal.size();
44 
45  // Initialise buddy check PGE to missing data indicator.
46  PGEBd.assign(obsVal.size(), PGEMDI);
47 
48  // Combined (obs and bkg) error variance.
49  float ErrVar = 0.0;
50  // Squared difference from background / ErrVar.
51  double SDiff = 0.0;
52  // Probability density of good observation.
53  double PdGood = 0.0;
54  // PGE after background check.
55  double PGEBk = 0.0;
56 
57  const bool obsErrEmpty = obsErr.empty();
58  const bool bkgErrEmpty = bkgErr.empty();
59 
60  for (size_t jloc = 0; jloc < numLocs; ++jloc) {
61  // Calculate combined error variance.
62  if (!obsErrEmpty && !bkgErrEmpty &&
63  obsErr[jloc] >= 0 && bkgErr[jloc] >= 0) {
64  ErrVar = std::pow(ObErrMult * obsErr[jloc], 2) +
65  std::pow(BkgErrMult * bkgErr[jloc], 2);
66  } else {
67  ErrVar = missingValueFloat;
68  }
69  // Set combined error variance to maximum value (if defined).
70  if (ErrVarMax > 0.0) {
71  ErrVar = std::min(ErrVar, ErrVarMax);
72  }
73 
74  // Update PGE (if all of the required values are present).
75  if (obsVal[jloc] != missingValueFloat &&
76  bkgVal[jloc] != missingValueFloat &&
77  ErrVar != missingValueFloat) {
78  if (obsVal2 && bkgVal2 &&
79  (*obsVal2)[jloc] != missingValueFloat &&
80  (*bkgVal2)[jloc] != missingValueFloat) { // Vector observable.
81  SDiff = (std::pow(obsVal[jloc] - bkgVal[jloc], 2) +
82  std::pow((*obsVal2)[jloc] - (*bkgVal2)[jloc], 2)) / static_cast<double>(ErrVar);
83  // Bivariate normal distribution; square root does not appear in denominator.
84  PdGood = std::exp(-0.5 * std::min(SDiff, 2.0 * ExpArgMax)) / (2.0 * M_PI * ErrVar);
85  } else { // Scalar observable.
86  SDiff = std::pow(obsVal[jloc] - bkgVal[jloc], 2) / ErrVar;
87  // Univariate normal distribution; square root appears in denominator.
88  PdGood = std::exp(-0.5 * std::min(SDiff, 2.0 * ExpArgMax)) /
89  std::sqrt(2.0 * M_PI * ErrVar);
90  }
91 
92  // Calculate PGE after background check, normalising appropriately.
93  PGEBk = (PdBad[jloc] * PGE[jloc]) /
94  (PdBad[jloc] * PGE[jloc] + PdGood * (1.0 - PGE[jloc]));
95 
96  // Set QC flags.
98  if (PGEBk >= PGECrit) {
100  }
101  } else {
102  // Deal with missing data.
103  SDiff = SDiffCrit;
104  PGEBk = PGEMDI;
105  }
106 
107  // Pack PGEs for use in later routines.
108  PGE[jloc] = trunc(PGEBk * PGEMult) + PGE[jloc];
109 
110  // Model-level data may have additional processing.
111  if (ModelLevels &&
112  (SDiff >= SDiffCrit ||
115  PGEBk = PGEMDI; // Do not apply buddy check in this case.
117  }
118 
119  // Update PGE for buddy check.
120  PGEBd[jloc] = PGEBk;
121  }
122  }
123 } // namespace ufo
Options controlling the operation of the calculations involving probability of gross error.
oops::Parameter< float > PGE_BkgErrMult
Multiplication factor for background errors.
oops::Parameter< float > PGE_PGECrit
PGE rejection limit.
oops::Parameter< float > PGE_SDiffCrit
Critical value for squared difference from background / ErrVar.
oops::Parameter< float > PGE_ExpArgMax
Maximum value of exponent in background QC.
oops::Parameter< float > PGE_ObErrMult
Multiplication factor for observation errors.
@ FinalRejectFlag
Final QC flag.
@ BackPerfFlag
Background check performed.
@ BackRejectFlag
PGE>0.5 after backgr check.
@ PermRejectFlag
Blacklisted data.
Definition: RunCRTM.h:27
void BayesianPGEUpdate(const ProbabilityOfGrossErrorParameters &options, const std::vector< float > &obsVal, const std::vector< float > &obsErr, const std::vector< float > &bkgVal, const std::vector< float > &bkgErr, const std::vector< float > &PdBad, const bool ModelLevels, std::vector< int > &flags, std::vector< float > &PGE, std::vector< float > &PGEBd, float ErrVarMax, const std::vector< float > *obsVal2, const std::vector< float > *bkgVal2)
Bayesian update of probability of gross error (PGE)