UFO
TrackCheckUtils.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 <algorithm>
9 #include <cmath>
10 #include <memory>
11 #include <string>
12 
13 #include "eckit/exception/Exceptions.h"
14 #include "eckit/geometry/Point2.h"
15 #include "eckit/geometry/Point3.h"
16 #include "eckit/geometry/Sphere.h"
17 #include "ioda/ObsDataVector.h"
18 #include "ioda/ObsSpace.h"
19 #include "oops/util/DateTime.h"
20 #include "oops/util/sqr.h"
21 #include "ufo/filters/QCflags.h"
24 #include "ufo/utils/Constants.h"
26 
27 namespace ufo {
28 namespace {
29 /// \brief Returns the square of the distance between the two \p Point arguments.
31  float sum = 0;
32  for (size_t i = 0; i < a.size(); ++i)
33  sum += util::sqr(a[i] - b[i]);
34  return sum;
35 }
36 
37 TrackCheckUtils::Point pointFromLatLon(float latitude, float longitude) {
38  // This local copy is needed because convertSphericalToCartesian takes the first parameter by
39  // reference, but Constants::mean_earth_rad has no out-of-line definition.
40  const double meanEarthRadius = Constants::mean_earth_rad;
41  eckit::geometry::Point3 eckitPoint;
42  eckit::geometry::Sphere::convertSphericalToCartesian(
43  meanEarthRadius, eckit::geometry::Point2(longitude, latitude), eckitPoint);
45  std::copy(eckitPoint.begin(), eckitPoint.end(), Point.begin());
46  return Point;
47 }
48 } // namespace
49 
50 /// \p Returns the distance between the two cartesian-mapped \p Point arguments
51 float TrackCheckUtils::distance(const Point &a, const Point &b) {
52  return std::sqrt(distance2(a, b));
53 }
54 
55 /// Return the vector of elements of \p categories with indices \p validObsIds.
56 template <typename T>
57 std::vector<T> getValidObservationCategories(const std::vector<T> &categories,
58  const std::vector<size_t> validObsIds) {
59  std::vector<T> validObsCategories(validObsIds.size());
60  for (size_t validObsIndex = 0; validObsIndex < validObsIds.size(); ++validObsIndex) {
61  validObsCategories[validObsIndex] = categories[validObsIds[validObsIndex]];
62  }
63  return validObsCategories;
64 }
65 
67  const std::vector<bool> &apply, const std::shared_ptr<ioda::ObsDataVector<int>> &flags) {
68  std::vector<size_t> validObsIds;
69  for (size_t obsId = 0; obsId < apply.size(); ++obsId)
70  if (apply[obsId] && (*(flags))[0][obsId] == QCflags::pass)
71  validObsIds.push_back(obsId);
72  return validObsIds;
73 }
74 
75 void TrackCheckUtils::groupObservationsByStation(const std::vector<size_t> &validObsIds,
76  RecursiveSplitter &splitter,
77  const eckit::Configuration &config,
78  const ioda::ObsSpace &obsdb) {
79  std::unique_ptr<TrackCheckUtilsParameters> baseOptions_;
80  baseOptions_.reset(new TrackCheckUtilsParameters());
81  baseOptions_->deserialize(config);
82  if (baseOptions_->stationIdVariable.value() == boost::none) {
83  if (obsdb.obs_group_var().empty()) {
84  // Observations were not grouped into records.
85  // Assume all observations were taken during the same station.
86  return;
87  } else {
88  groupObservationsByRecordNumber(validObsIds, splitter, obsdb);
89  }
90  } else {
91  groupObservationsByVariable(*baseOptions_->stationIdVariable.value(),
92  validObsIds, splitter, obsdb);
93  }
94 }
95 void TrackCheckUtils::groupObservationsByRecordNumber(const std::vector<size_t> &validObsIds,
96  RecursiveSplitter &splitter,
97  const ioda::ObsSpace &obsdb) {
98  const std::vector<size_t> &obsCategories = obsdb.recnum();
99  std::vector<size_t> validObsCategories = getValidObservationCategories(
100  obsCategories, validObsIds);
101  splitter.groupBy(validObsCategories);
102 }
103 
105  const std::vector<size_t> &validObsIds,
106  RecursiveSplitter &splitter,
107  const ioda::ObsSpace &obsdb) {
108  switch (obsdb.dtype(variable.group(), variable.variable())) {
109  case ioda::ObsDtype::Integer:
110  groupObservationsByTypedVariable<int>(variable, validObsIds, splitter, obsdb);
111  break;
112 
113  case ioda::ObsDtype::String:
114  groupObservationsByTypedVariable<std::string>(variable, validObsIds, splitter, obsdb);
115  break;
116 
117  default:
118  throw eckit::UserError(
119  "Only integer and string variables may be used as station IDs", Here());
120  }
121 }
122 template <typename VariableType>
124  const std::vector<size_t> &validObsIds,
125  RecursiveSplitter &splitter,
126  const ioda::ObsSpace &obsdb) {
127  std::vector<VariableType> obsCategories(obsdb.nlocs());
128  obsdb.get_db(variable.group(), variable.variable(), obsCategories);
129  std::vector<VariableType> validObsCategories = getValidObservationCategories(
130  obsCategories, validObsIds);
131 
132  splitter.groupBy(validObsCategories);
133 }
134 
135 void TrackCheckUtils::sortTracksChronologically(const std::vector<size_t> &validObsIds,
136  RecursiveSplitter &splitter,
137  const ioda::ObsSpace &obsdb) {
138  std::vector<util::DateTime> times(obsdb.nlocs());
139  obsdb.get_db("MetaData", "datetime", times);
140  splitter.sortGroupsBy([&times, &validObsIds](size_t obsIndexA, size_t obsIndexB)
141  { return times[validObsIds[obsIndexA]] < times[validObsIds[obsIndexB]]; });
142 }
143 
145  TrackCheckUtils::collectObservationsLocations(const ioda::ObsSpace &obsdb) {
146  ObsGroupLocationTimes locationTimes;
147 
148  locationTimes.latitudes.resize(obsdb.nlocs());
149  obsdb.get_db("MetaData", "latitude", locationTimes.latitudes);
150 
151  locationTimes.longitudes.resize(obsdb.nlocs());
152  obsdb.get_db("MetaData", "longitude", locationTimes.longitudes);
153 
154  locationTimes.datetimes.resize(obsdb.nlocs());
155  obsdb.get_db("MetaData", "datetime", locationTimes.datetimes);
156 
157  return locationTimes;
158 }
159 
160 void TrackCheckUtils::flagRejectedObservations(const std::vector<bool> &isRejected,
161  std::vector<std::vector<bool> >
162  &flagged) {
163  for (std::vector<bool> & variableFlagged : flagged)
164  for (size_t obsId = 0; obsId < isRejected.size(); ++obsId)
165  if (isRejected[obsId])
166  variableFlagged[obsId] = true;
167 }
168 
170  const util::DateTime &time)
171  : location_(pointFromLatLon(latitude, longitude)), time_(time)
172 {}
173 
175  return numChecks_ != 0 ? static_cast<float>(numFailedChecks_) / numChecks_ : 0.0f;
176 }
178  : numFailedChecks_(0), numChecks_(0)
179 {}
181  if (result != CheckResult::SKIPPED)
182  ++numChecks_;
183  if (result == CheckResult::FAILED)
184  ++numFailedChecks_;
185  assert(numChecks_ >= 0);
186  assert(numFailedChecks_ >= 0);
187  assert(numFailedChecks_ <= numChecks_);
188 }
189 
191  if (result != CheckResult::SKIPPED)
192  --numChecks_;
193  if (result == CheckResult::FAILED)
194  --numFailedChecks_;
195  assert(numChecks_ >= 0);
196  assert(numFailedChecks_ >= 0);
197  assert(numFailedChecks_ <= numChecks_);
198 }
199 
200 } // namespace ufo
ufo::anonymous_namespace{TrackCheckUtils.cc}::distance2
float distance2(const TrackCheckUtils::Point &a, const TrackCheckUtils::Point &b)
Returns the square of the distance between the two Point arguments.
Definition: TrackCheckUtils.cc:30
ufo::TrackCheckUtils::ObsLocationTime::ObsLocationTime
ObsLocationTime(float latitude, float longitude, const util::DateTime &time)
Definition: TrackCheckUtils.cc:169
ufo::RecursiveSplitter::sortGroupsBy
void sortGroupsBy(Compare comp)
Sort the elements in each equivalence class in ascending order.
Definition: src/ufo/utils/RecursiveSplitter.h:329
ufo::QCflags::pass
constexpr int pass
Definition: QCflags.h:14
ufo::TrackCheckUtils::CheckCounter::CheckCounter
CheckCounter()
Definition: TrackCheckUtils.cc:177
ufo::getValidObservationCategories
std::vector< T > getValidObservationCategories(const std::vector< T > &categories, const std::vector< size_t > validObsIds)
Return the vector of elements of categories with indices validObsIds.
Definition: TrackCheckUtils.cc:57
ufo::TrackCheckUtils::getValidObservationIds
std::vector< size_t > getValidObservationIds(const std::vector< bool > &apply, const std::shared_ptr< ioda::ObsDataVector< int >> &flags)
Definition: TrackCheckUtils.cc:66
ufo::TrackCheckUtils::groupObservationsByTypedVariable
void groupObservationsByTypedVariable(const Variable &variable, const std::vector< size_t > &validObsIds, RecursiveSplitter &splitter, const ioda::ObsSpace &obsdb)
Definition: TrackCheckUtils.cc:123
ufo::RecursiveSplitter
Partitions an array into groups of elements equivalent according to certain criteria.
Definition: src/ufo/utils/RecursiveSplitter.h:63
ufo::TrackCheckUtils::CheckCounter::failedChecksFraction
float failedChecksFraction() const
Definition: TrackCheckUtils.cc:174
ufo::TrackCheckUtils::CheckCounter::registerCheckResult
void registerCheckResult(const CheckResult &result)
Definition: TrackCheckUtils.cc:180
ufo::TrackCheckUtils::CheckResult::SKIPPED
@ SKIPPED
ufo::TrackCheckUtils::ObsGroupLocationTimes::datetimes
std::vector< util::DateTime > datetimes
Definition: TrackCheckUtils.h:55
ufo
Definition: RunCRTM.h:27
ufo::Constants::mean_earth_rad
static constexpr double mean_earth_rad
Definition: Constants.h:39
ufo::TrackCheckUtils::CheckCounter::unregisterCheckResult
void unregisterCheckResult(const CheckResult &result)
Definition: TrackCheckUtils.cc:190
ufo::TrackCheckUtils::groupObservationsByVariable
void groupObservationsByVariable(const Variable &variable, const std::vector< size_t > &validObsIds, RecursiveSplitter &splitter, const ioda::ObsSpace &obsdb)
Definition: TrackCheckUtils.cc:104
ufo::RecursiveSplitter::groupBy
void groupBy(const std::vector< size_t > &categories)
Split existing equivalence classes according to a new criterion.
Definition: src/ufo/utils/RecursiveSplitter.h:275
ufo::TrackCheckUtils::Point
std::array< float, 3 > Point
Definition: TrackCheckUtils.h:38
ufo::TrackCheckUtilsParameters
Options controlling the operation of the track check filter.
Definition: TrackCheckUtilsParameters.h:22
ufo::TrackCheckUtils::groupObservationsByStation
void groupObservationsByStation(const std::vector< size_t > &validObsIds, RecursiveSplitter &splitter, const eckit::Configuration &config, const ioda::ObsSpace &obsdb)
Definition: TrackCheckUtils.cc:75
TrackCheckUtilsParameters.h
ufo::anonymous_namespace{TrackCheckUtils.cc}::pointFromLatLon
TrackCheckUtils::Point pointFromLatLon(float latitude, float longitude)
Definition: TrackCheckUtils.cc:37
ufo::TrackCheckUtils::collectObservationsLocations
ObsGroupLocationTimes collectObservationsLocations(const ioda::ObsSpace &obsdb)
Definition: TrackCheckUtils.cc:145
QCflags.h
ufo::TrackCheckUtils::ObsGroupLocationTimes::latitudes
std::vector< float > latitudes
Definition: TrackCheckUtils.h:53
ufo::Variable::group
const std::string & group() const
Definition: Variable.cc:117
ufo::Variable::variable
const std::string & variable() const
Definition: Variable.cc:100
TrackCheckUtils.h
ioda::ObsDataVector< int >
ufo::TrackCheckUtils::ObsGroupLocationTimes
Locations/times of all observations processed by the track checking filter.
Definition: TrackCheckUtils.h:51
ufo::TrackCheckUtils::distance
float distance(const Point &a, const Point &b)
Returns the distance between the two cartesian-mapped Point arguments
Definition: TrackCheckUtils.cc:51
Constants.h
ufo::TrackCheckUtils::sortTracksChronologically
void sortTracksChronologically(const std::vector< size_t > &validObsIds, RecursiveSplitter &splitter, const ioda::ObsSpace &obsdb)
Definition: TrackCheckUtils.cc:135
ufo::TrackCheckUtils::groupObservationsByRecordNumber
void groupObservationsByRecordNumber(const std::vector< size_t > &validObsIds, RecursiveSplitter &splitter, const ioda::ObsSpace &obsdb)
Definition: TrackCheckUtils.cc:95
ufo::Variable
Definition: Variable.h:23
ufo::TrackCheckUtils::CheckResult::FAILED
@ FAILED
ufo::TrackCheckUtils::flagRejectedObservations
void flagRejectedObservations(const std::vector< bool > &isRejected, std::vector< std::vector< bool > > &flagged)
Definition: TrackCheckUtils.cc:160
RecursiveSplitter.h
ufo::TrackCheckUtils::ObsGroupLocationTimes::longitudes
std::vector< float > longitudes
Definition: TrackCheckUtils.h:54
ufo::TrackCheckUtils::CheckResult
CheckResult
Definition: TrackCheckUtils.h:42