UFO
MetOfficeBuddyCheck.cc
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2020 Met Office UK
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 
8 #include "ufo/filters/MetOfficeBuddyCheck.h"
9 
10 #include <algorithm>
11 #include <cmath>
12 #include <functional>
13 #include <iomanip>
14 #include <map>
15 #include <string>
16 #include <utility>
17 #include <vector>
18 
19 #include <boost/format.hpp>
20 
21 #include "eckit/config/Configuration.h"
22 #include "ioda/ObsDataVector.h"
23 #include "ioda/ObsSpace.h"
24 #include "oops/base/Variables.h"
25 #include "oops/util/DateTime.h"
26 #include "oops/util/Duration.h"
27 #include "oops/util/Logger.h"
28 #include "oops/util/sqr.h"
29 #include "ufo/filters/MetOfficeBuddyCheckParameters.h"
30 #include "ufo/filters/MetOfficeBuddyPair.h"
31 #include "ufo/filters/MetOfficeBuddyPairFinder.h"
32 #include "ufo/utils/PiecewiseLinearInterpolation.h"
33 
34 namespace ufo {
35 
36 namespace {
37 
38 const double expArgMax = 80.0; // maximum argument of the exponential function
39 const float maxGrossErrorProbability = 1.0; // maximum allowed value for PGE
40 
41 std::string fullVariableName(const Variable &var)
42 {
43  if (var.group().empty())
44  return var.variable();
45  else
46  return var.variable() + "@" + var.group();
47 }
48 
49 template <typename T>
50 std::vector<T> uniqueElements(const std::vector<T> &v)
51 {
52  std::vector<T> result(v);
53  std::sort(result.begin(), result.end());
54  const auto newEnd = std::unique(result.begin(), result.end());
55  result.erase(newEnd, result.end());
56  return result;
57 }
58 
59 template <typename T>
60 std::vector<int> mapDistinctValuesToDistinctInts(const std::vector<T> &values)
61 {
62  const std::vector<T> uniqueValues = uniqueElements(values);
63  std::map<T, int> intByValue;
64  {
65  int index = 0;
66  for (const T& value : uniqueValues)
67  intByValue[value] = index++;
68  }
69 
70  std::vector<int> ints;
71  ints.reserve(values.size());
72  std::transform(values.begin(), values.end(), std::back_inserter(ints),
73  [&](const T& value) { return intByValue.at(value); });
74  return ints;
75 }
76 
77 struct ScalarSingleLevelVariableData {
78  const std::vector<int> *varFlags = nullptr;
79  const std::vector<float> *obsValues = nullptr;
80  const std::vector<float> *obsBiases = nullptr;
81  const std::vector<float> *obsErrors = nullptr;
82  std::vector<float> bgValues;
83  std::vector<float> bgErrors;
84  std::vector<float> grossErrorProbabilities;
85 };
86 
87 } // namespace
88 
89 /// \brief Metadata of all observations processed by the filter.
90 struct MetOfficeBuddyCheck::MetaData {
91  std::vector<float> latitudes;
92  std::vector<float> longitudes;
93  std::vector<util::DateTime> datetimes;
94  boost::optional<std::vector<float>> pressures;
95  std::vector<int> stationIds;
96 };
97 
98 MetOfficeBuddyCheck::MetOfficeBuddyCheck(ioda::ObsSpace& obsdb, const eckit::Configuration& config,
99  std::shared_ptr<ioda::ObsDataVector<int> > flags,
100  std::shared_ptr<ioda::ObsDataVector<float> > obserr)
101  : FilterBase(obsdb, config, std::move(flags), std::move(obserr))
102 {
103  oops::Log::debug() << "MetOfficeBuddyCheck: config = " << config_ << std::endl;
104 
105  options_.reset(new MetOfficeBuddyCheckParameters());
106  options_->deserialize(config);
107 
108  allvars_ += Variables(filtervars_, "HofX");
109 }
110 
111 // Required for the correct destruction of options_.
112 MetOfficeBuddyCheck::~MetOfficeBuddyCheck()
113 {}
114 
115 void MetOfficeBuddyCheck::applyFilter(const std::vector<bool> & apply,
116  const Variables & filtervars,
117  std::vector<std::vector<bool>> & flagged) const {
118  // Identify observations to process and collect their data and metadata
119 
120  const std::vector<size_t> validObsIds = getValidObservationIds(apply);
121  MetaData obsData = collectMetaData();
122  const std::vector<float> bgErrorHorizCorrScales = calcBackgroundErrorHorizontalCorrelationScales(
123  validObsIds, obsData.latitudes);
124  const std::vector<bool> verbose = flagAndPrintVerboseObservations(
125  validObsIds, obsData.latitudes, obsData.longitudes, obsData.datetimes,
126  obsData.pressures.get_ptr(), obsData.stationIds, bgErrorHorizCorrScales);
127 
128  const oops::Variables &observedVars = obsdb_.obsvariables();
129  ioda::ObsDataVector<float> obsValues(obsdb_, filtervars.toOopsVariables(), "ObsValue");
130  ioda::ObsDataVector<float> obsBiases(obsdb_, filtervars.toOopsVariables(), "ObsBias",
131  false /*fail if not found?*/);
132 
133  auto getFilterVariableName = [&] (size_t filterVarIndex) {
134  return filtervars.variable(filterVarIndex).variable();
135  };
136 
137  auto getScalarSingleLevelVariableData = [&] (size_t filterVarIndex) {
138  const size_t observedVarIndex = observedVars.find(getFilterVariableName(filterVarIndex));
139 
140  ScalarSingleLevelVariableData data;
141  data.varFlags = &(*flags_)[observedVarIndex];
142  data.obsValues = &obsValues[filterVarIndex];
143  data.obsBiases = &obsBiases[filterVarIndex];
144  data.obsErrors = &(*obserr_)[observedVarIndex];
145  data_.get(ufo::Variable(filtervars[filterVarIndex], "HofX"), data.bgValues);
146  // TODO(wsmigaj): It isn't clear where to get background error variances from.
147  data_.get(ufo::Variable(filtervars[filterVarIndex], "HofXError"), data.bgErrors);
148  // TODO(wsmigaj): How is this variable going to be initialized?
149  data_.get(ufo::Variable(filtervars[filterVarIndex], "GrossErrorProbability"),
150  data.grossErrorProbabilities);
151  return data;
152  };
153 
154  // Identify buddy pairs
155 
156  const std::vector<float> *pressures =
157  options_->sortByPressure ? obsData.pressures.get_ptr() : nullptr;
158  MetOfficeBuddyPairFinder buddyPairFinder(*options_, obsData.latitudes, obsData.longitudes,
159  obsData.datetimes, pressures,
160  obsData.stationIds);
161  const std::vector<MetOfficeBuddyPair> buddyPairs = buddyPairFinder.findBuddyPairs(validObsIds);
162 
163  // Buddy-check all filter variables
164 
165  // Gross error probabilities updated by buddy check, indexed by variable name.
166  // (The updated values are copied to the ObsSpace only after all variables have been processed).
167  std::map<std::string, std::vector<float>> calculatedGrossErrProbsByVarName;
168 
169  bool previousVariableWasFirstComponentOfTwo = false;
170  for (size_t filterVarIndex = 0; filterVarIndex < filtervars.size(); ++filterVarIndex) {
171  if (previousVariableWasFirstComponentOfTwo) {
172  // Vector (two-component) variable
173  ScalarSingleLevelVariableData firstComponentData =
174  getScalarSingleLevelVariableData(filterVarIndex - 1);
175  ScalarSingleLevelVariableData secondComponentData =
176  getScalarSingleLevelVariableData(filterVarIndex);
177 
178  for (size_t i = 0; i < firstComponentData.grossErrorProbabilities.size(); ++i)
179  firstComponentData.grossErrorProbabilities[i] =
180  std::max(firstComponentData.grossErrorProbabilities[i],
181  secondComponentData.grossErrorProbabilities[i]);
182 
183  checkVectorSurfaceData(buddyPairs,
184  *firstComponentData.varFlags,
185  verbose, bgErrorHorizCorrScales,
186  obsData.stationIds, obsData.datetimes,
187  *firstComponentData.obsValues, *firstComponentData.obsBiases,
188  *secondComponentData.obsValues, *secondComponentData.obsBiases,
189  *firstComponentData.obsErrors,
190  firstComponentData.bgValues, secondComponentData.bgValues,
191  firstComponentData.bgErrors,
192  firstComponentData.grossErrorProbabilities);
193  // OPS doesn't update the gross error probabilities of the second component variable,
194  // but it seems more consistent to do so (and it facilitates the implementation of
195  // flagRejectedObservations().
196  // The implementation of checkVectorSurfaceData() still assumes that the *input* gross error
197  // probabilities, flags and background error estimates are the same for both components.
198  calculatedGrossErrProbsByVarName[getFilterVariableName(filterVarIndex - 1)] =
199  firstComponentData.grossErrorProbabilities;
200  calculatedGrossErrProbsByVarName[getFilterVariableName(filterVarIndex)] =
201  std::move(firstComponentData.grossErrorProbabilities);
202 
203  previousVariableWasFirstComponentOfTwo = false;
204  } else {
205  if (filtervars[filterVarIndex].options().getBool("first_component_of_two", false)) {
206  previousVariableWasFirstComponentOfTwo = true;
207  } else {
208  // Scalar variable
209  ScalarSingleLevelVariableData data =
210  getScalarSingleLevelVariableData(filterVarIndex);
211  checkScalarSurfaceData(buddyPairs,
212  *data.varFlags, verbose, bgErrorHorizCorrScales,
213  obsData.stationIds, obsData.datetimes,
214  *data.obsValues, *data.obsBiases, *data.obsErrors,
215  data.bgValues, data.bgErrors,
216  data.grossErrorProbabilities);
217  calculatedGrossErrProbsByVarName[getFilterVariableName(filterVarIndex)] =
218  std::move(data.grossErrorProbabilities);
219  }
220  }
221  }
222 
223  // Update observations flags and gross error probabilities
224 
225  for (const auto &varNameAndGrossErrProbs : calculatedGrossErrProbsByVarName) {
226  obsdb_.put_db("GrossErrorProbability", varNameAndGrossErrProbs.first,
227  varNameAndGrossErrProbs.second);
228  }
229 
230  flagRejectedObservations(filtervars, calculatedGrossErrProbsByVarName, flagged);
231 
232  if (filtervars.size() != 0) {
233  oops::Log::trace() << "MetOfficeBuddyCheck: flagged? = " << flagged[0] << std::endl;
234  }
235 }
236 
237 MetOfficeBuddyCheck::MetaData MetOfficeBuddyCheck::collectMetaData() const {
238  MetaData obsData;
239 
240  obsData.latitudes.resize(obsdb_.nlocs());
241  obsdb_.get_db("MetaData", "latitude", obsData.latitudes);
242 
243  obsData.longitudes.resize(obsdb_.nlocs());
244  obsdb_.get_db("MetaData", "longitude", obsData.longitudes);
245 
246  obsData.datetimes.resize(obsdb_.nlocs());
247  obsdb_.get_db("MetaData", "datetime", obsData.datetimes);
248 
249  if (obsdb_.has("MetaData", "air_pressure")) {
250  obsData.pressures = std::vector<float>(obsdb_.nlocs());
251  obsdb_.get_db("MetaData", "air_pressure", *obsData.pressures);
252  }
253 
254  obsData.stationIds = getStationIds();
255 
256  return obsData;
257 }
258 
259 std::vector<int> MetOfficeBuddyCheck::getStationIds() const {
260  const boost::optional<Variable> &stationIdVariable = options_->stationIdVariable.value();
261  if (stationIdVariable == boost::none) {
262  if (obsdb_.obs_group_var().empty()) {
263  // Observations were not grouped into records.
264  // Assume all observations were taken by the same station.
265  return std::vector<int>(obsdb_.nlocs(), 0);
266  } else {
267  const std::vector<size_t> &recordNumbers = obsdb_.recnum();
268  return std::vector<int>(recordNumbers.begin(), recordNumbers.end());
269  }
270  } else {
271  switch (obsdb_.dtype(stationIdVariable->group(), stationIdVariable->variable())) {
272  case ioda::ObsDtype::Integer:
273  {
274  std::vector<int> stationIds(obsdb_.nlocs());
275  obsdb_.get_db(stationIdVariable->group(), stationIdVariable->variable(), stationIds);
276  return stationIds;
277  }
278 
279  case ioda::ObsDtype::String:
280  {
281  std::vector<std::string> stringIds(obsdb_.nlocs());
282  obsdb_.get_db(stationIdVariable->group(), stationIdVariable->variable(), stringIds);
283  return mapDistinctValuesToDistinctInts(stringIds);
284  }
285 
286  default:
287  throw eckit::UserError("Only integer and string variables may be used as station IDs",
288  Here());
289  }
290  }
291 }
292 
293 std::vector<float> MetOfficeBuddyCheck::calcBackgroundErrorHorizontalCorrelationScales(
294  const std::vector<size_t> &validObsIds, const std::vector<float> &latitudes) const {
295 
296  std::vector<double> abscissas, ordinates;
297  for (const std::pair<const float, float> &xy :
298  options_->horizontalCorrelationScaleInterpolationPoints.value()) {
299  abscissas.push_back(xy.first);
300  ordinates.push_back(xy.second);
301  }
302  PiecewiseLinearInterpolation horizontalCorrelationScaleInterp(std::move(abscissas),
303  std::move(ordinates));
304 
305  std::vector<float> scales(latitudes.size());
306  for (size_t obsId : validObsIds) {
307  scales[obsId] = horizontalCorrelationScaleInterp(latitudes[obsId]);
308  }
309 
310  return scales;
311 }
312 
313 std::vector<bool> MetOfficeBuddyCheck::flagAndPrintVerboseObservations(
314  const std::vector<size_t> &validObsIds,
315  const std::vector<float> &latitudes,
316  const std::vector<float> &longitudes,
317  const std::vector<util::DateTime> &datetimes,
318  const std::vector<float> *pressures,
319  const std::vector<int> &stationIds,
320  const std::vector<float> &bgErrorHorizCorrScales) const {
321 
322  size_t numVerboseObs = 0;
323  std::vector<bool> verbose(latitudes.size(), false);
324 
325  for (size_t obsId : validObsIds) {
326  verbose[obsId] = std::any_of(
327  options_->tracedBoxes.value().begin(),
328  options_->tracedBoxes.value().end(),
329  [&](const LatLonBoxParameters &box)
330  { return box.contains(latitudes[obsId], longitudes[obsId]); });
331 
332  if (verbose[obsId]) {
333  if (numVerboseObs == 0) {
334  oops::Log::trace() << "Obs StationId Lat Lon Pressure "
335  "Datetime CorrScale\n";
336  }
337 
338  ++numVerboseObs;
339  const float pressure = pressures != nullptr ? (*pressures)[obsId] : 0;
340  oops::Log::trace() << boost::format("%5d %9d %7.2f %7.2f %8.0f %s %9.1f\n") %
341  obsId % stationIds[obsId] % latitudes[obsId] % longitudes[obsId] %
342  pressure % datetimes[obsId] % bgErrorHorizCorrScales[obsId];
343  }
344  }
345 
346  oops::Log::trace() << "Buddy check: " << numVerboseObs << " verbose observations" << std::endl;
347 
348  return verbose;
349 }
350 
351 void MetOfficeBuddyCheck::checkScalarSurfaceData(const std::vector<MetOfficeBuddyPair> &pairs,
352  const std::vector<int> &flags,
353  const std::vector<bool> &verbose,
354  const std::vector<float> &bgErrorHorizCorrScales,
355  const std::vector<int> &stationIds,
356  const std::vector<util::DateTime> &datetimes,
357  const std::vector<float> &obsValues,
358  const std::vector<float> &obsBiases,
359  const std::vector<float> &obsErrors,
360  const std::vector<float> &bgValues,
361  const std::vector<float> &bgErrors,
362  std::vector<float> &pges) const {
363  using util::sqr;
364 
365  const bool isMaster = obsdb_.comm().rank() == 0;
366  if (isMaster) {
367  oops::Log::trace() << __func__ << " "
368  << " dampingFactor1 = " << options_->dampingFactor1
369  << ", dampingFactor2 = " << options_->dampingFactor2 << '\n';
370  oops::Log::trace() << "ObsA ObsB StatIdA StatIdB DiffA DiffB "
371  "Dist Corr Agree PgeA PgeB Mult\n";
372  }
373 
374  const double invTemporalCorrScale = 1.0 / options_->temporalCorrelationScale.value().toSeconds();
375 
376  for (const MetOfficeBuddyPair &pair : pairs) {
377  const size_t jA = pair.obsIdA;
378  const size_t jB = pair.obsIdB;
379 
380  // Check that observations are valid and buddy check is required
381  if (!(flags[jA] == QCflags::pass && flags[jB] == QCflags::pass &&
382  pges[jA] < maxGrossErrorProbability && pges[jB] < maxGrossErrorProbability))
383  continue;
384 
385  // eqn 3.9
386  const double hcScale = 0.5 * (bgErrorHorizCorrScales[jA] + bgErrorHorizCorrScales[jB]);
387  const double scaledDist = pair.distanceInKm / hcScale;
388 
389  // Background error correlation between ob positions.
390  // Surface data; treat vertical correlation as 1.0
391  // eqns 3.10, 3.11
392  const double corr = (1.0 + scaledDist) *
393  std::exp(-scaledDist - sqr((datetimes[jA] - datetimes[jB]).toSeconds() *
394  invTemporalCorrScale));
395 
396  if (corr < 0.1)
397  continue; // skip to next pair
398 
399  // Differences from background
400  double diffA = obsValues[jA] + obsBiases[jA] - bgValues[jA];
401  double diffB = obsValues[jB] + obsBiases[jB] - bgValues[jB];
402  // Estimated error variances (ob+bk) (eqn 2.5)
403  double errVarA = sqr(obsErrors[jA]) + sqr(bgErrors[jA]);
404  double errVarB = sqr(obsErrors[jB]) + sqr(bgErrors[jB]);
405  // Background error covariance between ob positions (eqn 3.13)
406  double covar = corr * bgErrors[jA] * bgErrors[jB];
407  // (Total error correlation between ob positions)**2 (eqn 3.14)
408  double rho2 = sqr(covar) / (errVarA * errVarB);
409  // Argument for exponents
410  double expArg = -(0.5 * rho2 / (1.0 - rho2)) *
411  (sqr(diffA) / errVarA + sqr(diffB) / errVarB - 2.0 * diffA * diffB / covar);
412  expArg = options_->dampingFactor1 * (-0.5 * std::log(1.0 - rho2) + expArg); // exponent of
413  expArg = std::min(expArgMax, std::max(-expArgMax, expArg)); // eqn 3.18
414  // Z = P(OA)*P(OB)/P(OA and OB)
415  double z = 1.0 / (1.0 - (1.0 - pges[jA]) * (1.0 - pges[jB]) * (1.0 - std::exp(expArg)));
416  if (z <= 0.0)
417  z = 1.0; // rounding error control
418  z = std::pow(z, options_->dampingFactor2); // eqn 3.16
419  pges[jA] *= z; // eqn 3.17
420  pges[jB] *= z; // eqn 3.17
421  if (isMaster && (verbose[jA] || verbose[jB])) {
422  oops::Log::trace() << boost::format("%5d %5d %8d %8d "
423  "%5.1f %5.1f %6.1f "
424  "%5.3f %6.3f %6.3f %6.3f %6.3f\n") %
425  jA % jB % stationIds[jA] % stationIds[jB] %
426  diffA % diffB % pair.distanceInKm %
427  corr % std::exp(expArg) % pges[jA] % pges[jB] % z;
428  }
429  }
430 }
431 
432 void MetOfficeBuddyCheck::checkVectorSurfaceData(const std::vector<MetOfficeBuddyPair> &pairs,
433  const std::vector<int> &flags,
434  const std::vector<bool> &verbose,
435  const std::vector<float> &bgErrorHorizCorrScales,
436  const std::vector<int> &stationIds,
437  const std::vector<util::DateTime> &datetimes,
438  const std::vector<float> &uObsValues,
439  const std::vector<float> &uObsBiases,
440  const std::vector<float> &vObsValues,
441  const std::vector<float> &vObsBiases,
442  const std::vector<float> &obsErrors,
443  const std::vector<float> &uBgValues,
444  const std::vector<float> &vBgValues,
445  const std::vector<float> &bgErrors,
446  std::vector<float> &pges) const {
447  using util::sqr;
448 
449  const bool isMaster = obsdb_.comm().rank() == 0;
450  if (isMaster) {
451  oops::Log::trace() << __func__ << " "
452  << " dampingFactor1 = " << options_->dampingFactor1
453  << ", dampingFactor2 = " << options_->dampingFactor2 << '\n';
454  oops::Log::trace() << "ObsA ObsB StatIdA StatIdB LDiffA LDiffB TDiffA TDiffB "
455  "Dist Corr Agree PgeA PgeB Mult\n";
456  }
457 
458  const double invTemporalCorrScale = 1.0 / options_->temporalCorrelationScale.value().toSeconds();
459 
460  for (const MetOfficeBuddyPair &pair : pairs) {
461  const size_t jA = pair.obsIdA;
462  const size_t jB = pair.obsIdB;
463 
464  // Check that observations are valid and buddy check is required
465  if (!(flags[jA] == QCflags::pass && flags[jB] == QCflags::pass &&
466  pges[jA] < maxGrossErrorProbability && pges[jB] < maxGrossErrorProbability))
467  continue;
468 
469  // eqn 3.9
470  double horizCorrScale = 0.5 * (bgErrorHorizCorrScales[jA] + bgErrorHorizCorrScales[jB]);
471  double scaleDist = pair.distanceInKm / horizCorrScale;
472  // Background error correlation between ob positions.
473  // Surface data; treat vertical correlation as 1.0
474  // eqns 3.10, 3.11
475  const double lCorr = std::exp(-scaleDist - sqr((datetimes[jA] - datetimes[jB]).toSeconds() *
476  invTemporalCorrScale));
477 
478  if ((1.0 + scaleDist) * lCorr < 0.1)
479  continue; // skip to next pair
480 
481  // Calculate longitudinal and transverse wind components
482  double sinRot = std::sin(pair.rotationAInRad);
483  double cosRot = std::cos(pair.rotationAInRad);
484  // Difference from background - longitudinal wind
485  double lDiffA = cosRot * (uObsValues[jA] + uObsBiases[jA] - uBgValues[jA])
486  + sinRot * (vObsValues[jA] + vObsBiases[jA] - vBgValues[jA]); // eqn 3.19
487  // Difference from background - transverse wind
488  double tDiffA = - sinRot * (uObsValues[jA] + uObsBiases[jA] - uBgValues[jA])
489  + cosRot * (vObsValues[jA] + vObsBiases[jA] - vBgValues[jA]); // eqn 3.20
490  sinRot = std::sin(pair.rotationBInRad);
491  cosRot = std::cos(pair.rotationBInRad);
492  // Difference from background - longitudinal wind
493  double lDiffB = cosRot * (uObsValues[jB] + uObsBiases[jB] - uBgValues[jB])
494  + sinRot * (vObsValues[jB] + vObsBiases[jB] - vBgValues[jB]); // eqn 3.19
495  // Difference from background - transverse wind
496  double tDiffB = - sinRot * (uObsValues[jB] + uObsBiases[jB] - uBgValues[jB])
497  + cosRot * (vObsValues[jB] + vObsBiases[jB] - vBgValues[jB]); // eqn 3.20
498 
499  // Estimated error variances (ob + bk; component wind variance)
500  double errVarA = sqr(obsErrors[jA]) + sqr(bgErrors[jA]); // eqn 2.5
501  double errVarB = sqr(obsErrors[jB]) + sqr(bgErrors[jB]); // eqn 2.5
502 
503  // Calculate covariances and probabilities
504  double lCovar = lCorr * bgErrors[jA] * bgErrors[jB]; // eqn 3.13
505  double tCovar = (1.0 - options_->nonDivergenceConstraint * scaleDist) * lCovar; // eqn 3.12, 13
506  // rho2 = (total error correlation between ob positions)**2
507  double lRho2 = sqr(lCovar) / (errVarA * errVarB); // eqn 3.14
508  double tRho2 = sqr(tCovar) / (errVarA * errVarB); // eqn 3.14
509  // Argument for exponents
510  double expArg;
511  if (std::abs (tRho2) <= 0.00001)
512  expArg = 0.0; // prevent division by tCovar=0.0
513  else
514  expArg = -(0.5 * tRho2 / (1.0 - tRho2)) *
515  (sqr(tDiffA) / errVarA + sqr(tDiffB) / errVarB - 2.0 * tDiffA * tDiffB / tCovar);
516  expArg = expArg - (0.5 * lRho2 / (1.0 - lRho2)) *
517  (sqr(lDiffA) / errVarA + sqr(lDiffB) / errVarB - 2.0 * lDiffA * lDiffB / lCovar);
518  expArg = options_->dampingFactor1 * (-0.5 * std::log((1.0 - lRho2) * (1.0 - lRho2)) + expArg);
519  expArg = std::min(expArgMax, std::max(-expArgMax, expArg)); // eqn 3.22
520  // Z = P(OA)*P(OB)/P(OA and OB)
521  double z = 1.0 / (1.0 - (1.0 - pges[jA]) * (1.0 - pges[jB]) * (1.0 - std::exp(expArg)));
522  if (z <= 0.0)
523  z = 1.0; // rounding error control
524  z = std::pow(z, options_->dampingFactor2); // eqn 3.16
525  pges[jA] *= z; // eqn 3.17
526  pges[jB] *= z; // eqn 3.17
527 
528  if (isMaster && (verbose[jA] || verbose[jB])) {
529  oops::Log::trace() << boost::format("%5d %5d %8d %8d "
530  "%6.1f %6.1f %6.1f %6.1f %6.1f "
531  "%5.3f %6.3f %6.3f %6.3f %6.3f\n") %
532  jA % jB % stationIds[jA] % stationIds[jB] %
533  lDiffA % lDiffB % tDiffA % tDiffB % pair.distanceInKm %
534  lCorr % std::exp(expArg) % pges[jA] % pges[jB] % z;
535  }
536  }
537 }
538 
539 std::vector<size_t> MetOfficeBuddyCheck::getValidObservationIds(
540  const std::vector<bool> & apply) const {
541  std::vector<size_t> validObsIds;
542  for (size_t obsId = 0; obsId < apply.size(); ++obsId)
543  // TODO(wsmigaj): The second condition below may need reviewing.
544  // Perhaps we should process only observations marked as passed in all filter variables?
545  // Or those marked as passed in at least one filter variable?
546  if (apply[obsId] && (*flags_)[0][obsId] == QCflags::pass)
547  validObsIds.push_back(obsId);
548  return validObsIds;
549 }
550 
551 void MetOfficeBuddyCheck::flagRejectedObservations(
552  const Variables &filtervars,
553  const std::map<std::string, std::vector<float>> &grossErrProbsByVarName,
554  std::vector<std::vector<bool>> &flagged) const {
555  ASSERT(filtervars.size() == flagged.size());
556 
557  for (size_t varIndex = 0; varIndex < filtervars.size(); ++varIndex) {
558  const std::string varName = fullVariableName(filtervars[varIndex]);
559  const std::vector<float> &grossErrProbs = grossErrProbsByVarName.at(varName);
560  std::vector<bool> &variableFlagged = flagged[varIndex];
561  ASSERT(grossErrProbs.size() == variableFlagged.size());
562 
563  for (size_t obsId = 0; obsId < grossErrProbs.size(); ++obsId)
564  if (grossErrProbs[obsId] >= options_->rejectionThreshold)
565  variableFlagged[obsId] = true;
566  }
567 }
568 
569 void MetOfficeBuddyCheck::print(std::ostream & os) const {
570  os << "MetOfficeBuddyCheck: config = " << config_ << std::endl;
571 }
572 
573 } // namespace ufo
ufo::MetOfficeBuddyCheck::print
void print(std::ostream &) const override
Definition: MetOfficeBuddyCheck.cc:569
ufo::MetOfficeBuddyCheck::MetaData
Metadata of all observations processed by the filter.
Definition: MetOfficeBuddyCheck.cc:90
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::expArgMax
const double expArgMax
Definition: MetOfficeBuddyCheck.cc:38
MetOfficeBuddyCheck.h
ufo::QCflags::pass
constexpr int pass
Definition: QCflags.h:14
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::fullVariableName
std::string fullVariableName(const Variable &var)
Definition: MetOfficeBuddyCheck.cc:41
ufo::Variables
Definition: src/ufo/filters/Variables.h:24
MetOfficeBuddyPair.h
ufo::MetOfficeBuddyCheck::calcBackgroundErrorHorizontalCorrelationScales
std::vector< float > calcBackgroundErrorHorizontalCorrelationScales(const std::vector< size_t > &validObsIds, const std::vector< float > &latitudes) const
Calculates and returns background error correlation scales at observation locations.
Definition: MetOfficeBuddyCheck.cc:293
PiecewiseLinearInterpolation.h
ufo::MetOfficeBuddyCheck::options_
std::unique_ptr< MetOfficeBuddyCheckParameters > options_
Definition: src/ufo/filters/MetOfficeBuddyCheck.h:217
ufo_radiancerttov_utils_mod::debug
logical, public debug
Definition: ufo_radiancerttov_utils_mod.F90:100
ufo::FilterBase::obsdb_
ioda::ObsSpace & obsdb_
Definition: FilterBase.h:59
ufo::LatLonBoxParameters
A box covering a specified (closed) interval of latitudes and longitudes.
Definition: MetOfficeBuddyCheckParameters.h:25
ufo::MetOfficeBuddyCheck::checkScalarSurfaceData
void checkScalarSurfaceData(const std::vector< MetOfficeBuddyPair > &pairs, const std::vector< int > &flags, const std::vector< bool > &verbose, const std::vector< float > &bgErrorHorizCorrScales, const std::vector< int > &stationIds, const std::vector< util::DateTime > &datetimes, const std::vector< float > &obsValues, const std::vector< float > &obsBiases, const std::vector< float > &obsErrors, const std::vector< float > &bgValues, const std::vector< float > &bgErrors, std::vector< float > &pges) const
Buddy check for scalar surface quantities.
Definition: MetOfficeBuddyCheck.cc:351
ufo::FilterBase
FilterBase: Base class for UFO QC filters.
Definition: FilterBase.h:42
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::uniqueElements
std::vector< T > uniqueElements(const std::vector< T > &v)
Definition: MetOfficeBuddyCheck.cc:50
ufo::MetOfficeBuddyCheck::getValidObservationIds
std::vector< size_t > getValidObservationIds(const std::vector< bool > &apply) const
Returns a vector of IDs of all observations that should be buddy-checked.
Definition: MetOfficeBuddyCheck.cc:539
ufo::MetOfficeBuddyCheck::MetaData::latitudes
std::vector< float > latitudes
Definition: MetOfficeBuddyCheck.cc:91
ufo
Definition: RunCRTM.h:27
ufo::MetOfficeBuddyCheck::getStationIds
std::vector< int > getStationIds() const
Returns a vector of integer-valued station IDs, obtained from the source indicated by the filter para...
Definition: MetOfficeBuddyCheck.cc:259
MetOfficeBuddyPairFinder.h
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::ScalarSingleLevelVariableData::grossErrorProbabilities
std::vector< float > grossErrorProbabilities
Definition: MetOfficeBuddyCheck.cc:84
ufo::MetOfficeBuddyPair
Properties of a pair of observations checked against each other during buddy check.
Definition: MetOfficeBuddyPair.h:14
ufo::FilterBase::data_
ObsFilterData data_
Definition: FilterBase.h:65
ufo::MetOfficeBuddyCheckParameters
Options controlling the operation of the MetOfficeBuddyCheck filter.
Definition: MetOfficeBuddyCheckParameters.h:41
ufo::FilterBase::allvars_
ufo::Variables allvars_
Definition: FilterBase.h:63
ufo::Variables::size
size_t size() const
Definition: Variables.cc:92
ufo::PiecewiseLinearInterpolation
Represents a piecewise linear interpolation of a set of data points.
Definition: src/ufo/utils/PiecewiseLinearInterpolation.h:17
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::ScalarSingleLevelVariableData
Definition: MetOfficeBuddyCheck.cc:77
ufo::MetOfficeBuddyPairFinder
Finds pairs of close observations ("buddies") to check against each other.
Definition: src/ufo/filters/MetOfficeBuddyPairFinder.h:26
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::ScalarSingleLevelVariableData::bgValues
std::vector< float > bgValues
Definition: MetOfficeBuddyCheck.cc:82
ufo::FilterBase::flags_
std::shared_ptr< ioda::ObsDataVector< int > > flags_
Definition: FilterBase.h:61
ufo::MetOfficeBuddyCheck::flagAndPrintVerboseObservations
std::vector< bool > flagAndPrintVerboseObservations(const std::vector< size_t > &validObsIds, const std::vector< float > &latitudes, const std::vector< float > &longitudes, const std::vector< util::DateTime > &times, const std::vector< float > *pressures, const std::vector< int > &stationIds, const std::vector< float > &bgErrorHorizCorrScales) const
Identifies observations whose buddy checks should be logged.
Definition: MetOfficeBuddyCheck.cc:313
ufo::MetOfficeBuddyCheck::applyFilter
void applyFilter(const std::vector< bool > &, const Variables &, std::vector< std::vector< bool >> &) const override
Definition: MetOfficeBuddyCheck.cc:115
ufo::MetOfficeBuddyCheck::MetaData::pressures
boost::optional< std::vector< float > > pressures
Definition: MetOfficeBuddyCheck.cc:94
ufo::Variables::toOopsVariables
oops::Variables toOopsVariables() const
Definition: Variables.cc:144
ufo::Variable::group
const std::string & group() const
Definition: Variable.cc:117
ufo::abs
util::Duration abs(const util::Duration &duration)
Definition: TrackCheckUtils.h:31
ufo::Variable::variable
const std::string & variable() const
Definition: Variable.cc:100
ufo::MetOfficeBuddyCheck::MetaData::stationIds
std::vector< int > stationIds
Definition: MetOfficeBuddyCheck.cc:95
ufo::MetOfficeBuddyCheck::collectMetaData
MetaData collectMetaData() const
Collects and return smetadata of all observations.
Definition: MetOfficeBuddyCheck.cc:237
ufo::FilterBase::filtervars_
ufo::Variables filtervars_
Definition: FilterBase.h:64
ufo::MetOfficeBuddyCheck::flagRejectedObservations
void flagRejectedObservations(const Variables &filtervars, const std::map< std::string, std::vector< float >> &grossErrProbsByVarName, std::vector< std::vector< bool >> &flagged) const
Definition: MetOfficeBuddyCheck.cc:551
ufo::MetOfficeBuddyCheck::~MetOfficeBuddyCheck
~MetOfficeBuddyCheck() override
Definition: MetOfficeBuddyCheck.cc:112
ioda::ObsDataVector< int >
ufo::MetOfficeBuddyCheck::MetOfficeBuddyCheck
MetOfficeBuddyCheck(ioda::ObsSpace &obsdb, const eckit::Configuration &config, std::shared_ptr< ioda::ObsDataVector< int > > flags, std::shared_ptr< ioda::ObsDataVector< float > > obserr)
Definition: MetOfficeBuddyCheck.cc:98
ufo::Variables::variable
Variable variable(const size_t) const
Definition: Variables.cc:114
ufo::FilterBase::config_
const eckit::LocalConfiguration config_
Definition: FilterBase.h:60
ufo::MetOfficeBuddyCheck::MetaData::datetimes
std::vector< util::DateTime > datetimes
Definition: MetOfficeBuddyCheck.cc:93
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::ScalarSingleLevelVariableData::bgErrors
std::vector< float > bgErrors
Definition: MetOfficeBuddyCheck.cc:83
MetOfficeBuddyCheckParameters.h
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::mapDistinctValuesToDistinctInts
std::vector< int > mapDistinctValuesToDistinctInts(const std::vector< T > &values)
Definition: MetOfficeBuddyCheck.cc:60
ufo::ObsFilterData::get
void get(const Variable &, std::vector< float > &) const
Gets requested data from ObsFilterData.
Definition: ObsFilterData.cc:130
ufo::MetOfficeBuddyPairFinder::findBuddyPairs
std::vector< MetOfficeBuddyPair > findBuddyPairs(const std::vector< size_t > &validObsIds)
Returns a list of MetOfficeBuddyPair objects representing pairs of "buddy" observations that should b...
Definition: MetOfficeBuddyPairFinder.cc:40
ufo::Variable
Definition: Variable.h:23
ufo::MetOfficeBuddyCheck::checkVectorSurfaceData
void checkVectorSurfaceData(const std::vector< MetOfficeBuddyPair > &pairs, const std::vector< int > &flags, const std::vector< bool > &verbose, const std::vector< float > &bgErrorHorizCorrScales, const std::vector< int > &stationIds, const std::vector< util::DateTime > &datetimes, const std::vector< float > &uObsValues, const std::vector< float > &uObsBiases, const std::vector< float > &vObsValues, const std::vector< float > &vObsBiases, const std::vector< float > &obsErrors, const std::vector< float > &uBgValues, const std::vector< float > &vBgValues, const std::vector< float > &bgErrors, std::vector< float > &pges) const
Buddy check for vector (two-dimensional) surface quantities.
Definition: MetOfficeBuddyCheck.cc:432
ufo::MetOfficeBuddyCheck::MetaData::longitudes
std::vector< float > longitudes
Definition: MetOfficeBuddyCheck.cc:92
ufo::anonymous_namespace{MetOfficeBuddyCheck.cc}::maxGrossErrorProbability
const float maxGrossErrorProbability
Definition: MetOfficeBuddyCheck.cc:39