UFO
test/ufo/MetOfficeBuddyPairFinder.h
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 #ifndef TEST_UFO_METOFFICEBUDDYPAIRFINDER_H_
9 #define TEST_UFO_METOFFICEBUDDYPAIRFINDER_H_
10 
11 #include <iomanip>
12 #include <map>
13 #include <set>
14 #include <string>
15 #include <utility>
16 #include <vector>
17 
18 #include <boost/optional.hpp>
19 
20 #include "eckit/config/LocalConfiguration.h"
21 #include "eckit/testing/Test.h"
22 #include "ioda/ObsSpace.h"
23 #include "ioda/ObsVector.h"
24 #include "oops/mpi/mpi.h"
25 #include "oops/runs/Test.h"
26 #include "oops/util/Expect.h"
27 #include "test/TestEnvironment.h"
28 #include "ufo/filters/MetOfficeBuddyCheckParameters.h"
29 #include "ufo/filters/MetOfficeBuddyPairFinder.h"
30 #include "ufo/utils/StringUtils.h"
31 
32 namespace ufo {
33 namespace test {
34 
35 typedef std::pair<int, int> ObsPair;
36 
37 bool duplicatePairsPresent(const std::vector<ObsPair> &pairs) {
38  bool duplicatesFound = false;
39  std::set<ObsPair> pairSet;
40 
41  // Test for exact duplicates
42  for (const ObsPair &pair : pairs) {
43  bool inserted = pairSet.insert(pair).second;
44  if (!inserted) {
45  duplicatesFound = true;
46  oops::Log::trace() << "Duplicate pair found: (" << pair.first << ", " << pair.second << ")\n";
47  }
48  }
49 
50  // Test for duplicates in reverse order
51  for (const ObsPair &pair : pairs) {
52  ObsPair reversePair(pair.second, pair.first);
53  bool inserted = pairSet.insert(reversePair).second;
54  if (!inserted) {
55  duplicatesFound = true;
56  oops::Log::trace() << "Pair (" << pair.first << ", " << pair.second
57  << ") present also in reverse order\n";
58  }
59  }
60 
61  oops::Log::trace().flush();
62 
63  return duplicatesFound;
64 }
65 
66 template <typename Key, typename Value>
67 Value maxValue(const std::map<Key, Value> &map)
68 {
69  if (map.empty())
70  return Value();
71 
72  typedef typename std::map<Key, Value>::value_type ValueType;
73  return std::max_element(map.begin(), map.end(),
74  [](const ValueType &a, const ValueType &b)
75  { return a.second < b.second; })->second;
76 }
77 
78 int maxTotalNumBuddies(const std::vector<ObsPair> &pairs) {
79  std::map<int, int> numBuddiesByObsId;
80  for (const ObsPair & pair : pairs)
81  ++numBuddiesByObsId[pair.first];
82  return maxValue(numBuddiesByObsId);
83 }
84 
85 int maxNumBuddiesWithSameStationId(const std::vector<ObsPair> &pairs,
86  const std::vector<int> &stationIds) {
87  std::map<int, int> numBuddiesWithSameStationIdByObsId;
88  for (const ObsPair & pair : pairs)
89  if (stationIds[pair.first] == stationIds[pair.second])
90  ++numBuddiesWithSameStationIdByObsId[pair.first];
91  return maxValue(numBuddiesWithSameStationIdByObsId);
92 }
93 
94 void testDuplicatesAndBuddyCountConstraints(const eckit::LocalConfiguration &conf) {
95  util::DateTime bgn(conf.getString("window begin"));
96  util::DateTime end(conf.getString("window end"));
97 
98  const eckit::LocalConfiguration obsSpaceConf(conf, "obs space");
99  ioda::ObsSpace obsSpace(obsSpaceConf, oops::mpi::world(), bgn, end, oops::mpi::myself());
100 
101  boost::optional<std::vector<float>> airPressures;
102  if (obsSpace.has("MetaData", "air_pressure")) {
103  airPressures = std::vector<float>(obsSpace.nlocs());
104  obsSpace.get_db("MetaData", "air_pressure", *airPressures);
105  }
106 
107  std::vector<float> latitudes(obsSpace.nlocs());
108  obsSpace.get_db("MetaData", "latitude", latitudes);
109 
110  std::vector<float> longitudes(obsSpace.nlocs());
111  obsSpace.get_db("MetaData", "longitude", longitudes);
112 
113  std::vector<util::DateTime> datetimes(obsSpace.nlocs());
114  obsSpace.get_db("MetaData", "datetime", datetimes);
115 
116  std::vector<int> stationIds(obsSpace.recnum().begin(), obsSpace.recnum().end());
117 
118  std::vector<size_t> validObsIds;
119  if (conf.has("valid_obs_ids")) {
120  validObsIds = conf.getUnsignedVector("valid_obs_ids");
121  } else {
122  validObsIds.resize(obsSpace.nlocs());
123  std::iota(validObsIds.begin(), validObsIds.end(), 0);
124  }
125 
126  const eckit::LocalConfiguration modernFilterConf(conf, "Met Office Buddy Check, modern");
127  MetOfficeBuddyCheckParameters modernOptions;
128  modernOptions.deserialize(modernFilterConf);
129 
130  std::vector<ObsPair> modernPairs;
131  {
132  MetOfficeBuddyPairFinder finder(modernOptions, latitudes, longitudes, datetimes,
133  airPressures.get_ptr(), stationIds);
134  const std::vector<MetOfficeBuddyPair> buddyPairs = finder.findBuddyPairs(validObsIds);
135 
136  for (const MetOfficeBuddyPair & pair : buddyPairs) {
137  modernPairs.push_back(ObsPair(pair.obsIdA, pair.obsIdB));
138  }
139  }
140 
141  EXPECT(maxTotalNumBuddies(modernPairs) <= modernOptions.maxTotalNumBuddies);
142  EXPECT(maxNumBuddiesWithSameStationId(modernPairs, stationIds) <=
143  modernOptions.maxNumBuddiesWithSameStationId);
144  EXPECT_NOT(duplicatePairsPresent(modernPairs));
145 
146  const eckit::LocalConfiguration legacyFilterConf(conf, "Met Office Buddy Check, legacy");
147  MetOfficeBuddyCheckParameters legacyOptions;
148  legacyOptions.deserialize(legacyFilterConf);
149 
150  std::vector<ObsPair> legacyPairs;
151  {
152  MetOfficeBuddyPairFinder finder(legacyOptions, latitudes, longitudes, datetimes,
153  airPressures.get_ptr(), stationIds);
154  const std::vector<MetOfficeBuddyPair> buddyPairs = finder.findBuddyPairs(validObsIds);
155 
156  for (const MetOfficeBuddyPair & pair : buddyPairs) {
157  legacyPairs.push_back(ObsPair(pair.obsIdA, pair.obsIdB));
158  }
159  }
160 
161  EXPECT(maxTotalNumBuddies(legacyPairs) <= legacyOptions.maxTotalNumBuddies);
162  EXPECT(maxNumBuddiesWithSameStationId(legacyPairs, stationIds) <=
163  legacyOptions.maxNumBuddiesWithSameStationId);
164  EXPECT_NOT(duplicatePairsPresent(legacyPairs));
165 
166  // The legacy buddy collector sometimes fails to find all observations that can be classified
167  // as buddies while respecting constraints on their maximum number.
168  EXPECT(modernPairs.size() > legacyPairs.size());
169 }
170 
171 CASE("ufo/MetOfficeBuddyPairFinder/"
172  "Duplicates, constraints on buddy counts, legacy pair collector") {
173  testDuplicatesAndBuddyCountConstraints(eckit::LocalConfiguration(
174  ::test::TestEnvironment::config(),
175  "Duplicates, constraints on buddy counts, "
176  "legacy pair collector"));
177 }
178 
179 std::vector<MetOfficeBuddyPair> findBuddyPairs(const MetOfficeBuddyCheckParameters &options,
180  const std::vector<float> &latitudes,
181  const std::vector<float> &longitudes,
182  const std::vector<util::DateTime> &datetimes,
183  const std::vector<float> *pressures,
184  const std::vector<int> &stationIds,
185  const std::vector<size_t> &validObsIds) {
186  MetOfficeBuddyPairFinder finder(options, latitudes, longitudes, datetimes,
187  pressures, stationIds);
188  return finder.findBuddyPairs(validObsIds);
189 }
190 
191 void testInvarianceToLongitude(const eckit::LocalConfiguration &conf) {
192  const float searchRadius = 100; // km
193 
194  std::vector<float> referenceLatitudes = conf.getFloatVector("reference.latitudes");
195  std::vector<float> referenceLongitudes = conf.getFloatVector("reference.longitudes");
196 
197  std::vector<util::DateTime> datetimes(2, util::DateTime());
198  std::vector<int> stationIds{1, 2};
199 
200  std::vector<size_t> validObsIds{0, 1};
201 
202  const eckit::LocalConfiguration filterConf(conf, "Met Office Buddy Check");
203  MetOfficeBuddyCheckParameters options;
204  options.deserialize(filterConf);
205 
206  const std::vector<MetOfficeBuddyPair> referenceBuddyPairs = findBuddyPairs(
207  options, referenceLatitudes, referenceLongitudes,
208  datetimes, nullptr, stationIds, validObsIds);
209  EXPECT_EQUAL(referenceBuddyPairs.size(), 1);
210  const MetOfficeBuddyPair &refPair = referenceBuddyPairs[0];
211 
212  for (const eckit::LocalConfiguration &testConf : conf.getSubConfigurations("test")) {
213  std::vector<float> latitudes = testConf.getFloatVector("latitudes");
214  std::vector<float> longitudes = testConf.getFloatVector("longitudes");
215  const std::vector<MetOfficeBuddyPair> buddyPairs = findBuddyPairs(
216  options, latitudes, longitudes, datetimes, nullptr, stationIds, validObsIds);
217  EXPECT_EQUAL(buddyPairs.size(), 1);
218  const MetOfficeBuddyPair &pair = buddyPairs[0];
219 
220  if (pair.obsIdA == refPair.obsIdA) {
221  EXPECT_EQUAL(pair.obsIdB, refPair.obsIdB);
222  EXPECT(oops::is_close_relative(pair.distanceInKm, refPair.distanceInKm, 1e-6));
223  if (latitudes[pair.obsIdB] - latitudes[pair.obsIdA] ==
224  referenceLatitudes[refPair.obsIdB] - referenceLatitudes[refPair.obsIdA] &&
225  longitudes[pair.obsIdB] - longitudes[pair.obsIdA] ==
226  referenceLongitudes[refPair.obsIdB] - referenceLongitudes[refPair.obsIdA]) {
227  EXPECT(oops::is_close_relative(pair.rotationAInRad, refPair.rotationAInRad, 1e-6));
228  EXPECT(oops::is_close_relative(pair.rotationBInRad, refPair.rotationBInRad, 1e-6));
229  }
230  } else {
231  EXPECT_EQUAL(pair.obsIdA, refPair.obsIdB);
232  EXPECT_EQUAL(pair.obsIdB, refPair.obsIdA);
233  EXPECT(oops::is_close_relative(pair.distanceInKm, refPair.distanceInKm, 1e-6));
234  if (latitudes[pair.obsIdB] - latitudes[pair.obsIdA] ==
235  referenceLatitudes[refPair.obsIdA] - referenceLatitudes[refPair.obsIdB] &&
236  longitudes[pair.obsIdB] - longitudes[pair.obsIdA] ==
237  referenceLongitudes[refPair.obsIdA] - referenceLongitudes[refPair.obsIdB]) {
238  EXPECT(oops::is_close_relative(pair.rotationAInRad, refPair.rotationBInRad, 1e-6));
239  EXPECT(oops::is_close_relative(pair.rotationBInRad, refPair.rotationAInRad, 1e-6));
240  }
241  }
242  }
243 }
244 
245 CASE("ufo/TemporalThinning/Invariance to longitude, different zonal bands") {
246  testInvarianceToLongitude(eckit::LocalConfiguration(
247  ::test::TestEnvironment::config(),
248  "Invariance to longitude, different zonal bands"));
249 }
250 
251 CASE("ufo/TemporalThinning/Invariance to longitude, same zonal band") {
252  testInvarianceToLongitude(eckit::LocalConfiguration(
253  ::test::TestEnvironment::config(),
254  "Invariance to longitude, same zonal band"));
255 }
256 
257 void findEndpoint(double startLat, double startLon, double distance, double bearing,
258  double &endLat, double &endLon)
259 {
260  startLat *= ufo::Constants::deg2rad;
261  startLon *= ufo::Constants::deg2rad;
262  bearing *= ufo::Constants::deg2rad;
263  distance /= ufo::Constants::mean_earth_rad;
264  endLat = std::asin(std::sin(startLat) * std::cos(distance) +
265  std::cos(startLat) * std::sin(distance) * std::cos(bearing));
266  endLon = startLon + std::atan2(std::sin(bearing) * std::sin(distance) * std::cos(startLat),
267  std::cos(distance) - std::sin(startLat) * std::sin(endLat));
268  endLat *= ufo::Constants::rad2deg;
269  endLon *= ufo::Constants::rad2deg;
270 }
271 
272 void findEndpoint(float startLat, float startLon, float distance, float bearing,
273  float &endLat, float &endLon)
274 {
275  double dEndLat, dEndLon;
276  findEndpoint(startLat, startLon, distance, bearing, dEndLat, dEndLon);
277  endLon = dEndLon;
278  endLat = dEndLat;
279 }
280 
281 template <typename T>
282 bool contains(const std::set<T> &set, const T &element) {
283  return set.find(element) != set.end();
284 }
285 
286 ObsPair reverse(const ObsPair &pair) {
287  return ObsPair(pair.second, pair.first);
288 }
289 
290 void testSearchRadius(const eckit::LocalConfiguration &conf) {
291  const eckit::LocalConfiguration filterConf(conf, "Met Office Buddy Check");
292  MetOfficeBuddyCheckParameters options;
293  options.deserialize(filterConf);
294 
295  const float searchRadius = options.searchRadius;
296 
297  for (const eckit::LocalConfiguration &testConf : conf.getSubConfigurations("test")) {
298  float centralLatitude = testConf.getFloat("central_latitude");
299  float centralLongitude = testConf.getFloat("central_longitude");
300 
301  std::vector<float> latitudes{centralLatitude};
302  std::vector<float> longitudes{centralLongitude};
303  std::vector<ObsPair> pairsExpectedToBePresent, pairsExpectedToBeAbsent;
304 
305  const int numBearings = 8;
306  for (float i = 0; i < numBearings; ++i) {
307  float latitude, longitude;
308 
309  findEndpoint(centralLatitude, centralLongitude,
310  0.99f * searchRadius, i * 360.0f / numBearings, latitude, longitude);
311  latitudes.push_back(latitude);
312  longitudes.push_back(longitude);
313  pairsExpectedToBePresent.push_back(ObsPair(0, latitudes.size() - 1));
314 
315  findEndpoint(centralLatitude, centralLongitude,
316  1.01f * searchRadius, i * 360.0f / numBearings, latitude, longitude);
317  latitudes.push_back(latitude);
318  longitudes.push_back(longitude);
319  pairsExpectedToBeAbsent.push_back(ObsPair(0, latitudes.size() - 1));
320  }
321 
322  const std::vector<util::DateTime> datetimes(latitudes.size(), util::DateTime());
323  const std::vector<int> stationIds(latitudes.size(), 0);
324 
325  std::vector<size_t> validObsIds(latitudes.size());
326  std::iota(validObsIds.begin(), validObsIds.end(), 0);
327 
328  const std::vector<MetOfficeBuddyPair> buddyPairs = findBuddyPairs(
329  options, latitudes, longitudes, datetimes, nullptr, stationIds, validObsIds);
330 
331  std::set<ObsPair> pairs;
332  for (const MetOfficeBuddyPair & pair : buddyPairs) {
333  pairs.insert(ObsPair(pair.obsIdA, pair.obsIdB));
334  }
335 
336  for (const ObsPair &pair : pairsExpectedToBePresent)
337  EXPECT(contains(pairs, pair) || contains(pairs, reverse(pair)));
338 
339  for (const ObsPair &pair : pairsExpectedToBeAbsent)
340  EXPECT(!contains(pairs, pair) && !contains(pairs, reverse(pair)));
341  }
342 }
343 
344 CASE("ufo/TemporalThinning/Search radius") {
345  testSearchRadius(eckit::LocalConfiguration(::test::TestEnvironment::config(), "Search radius"));
346 }
347 
348 class MetOfficeBuddyPairFinder : public oops::Test {
349  private:
350  std::string testid() const override {return "ufo::test::MetOfficeBuddyPairFinder";}
351 
352  void register_tests() const override {}
353 
354  void clear() const override {}
355 };
356 
357 } // namespace test
358 } // namespace ufo
359 
360 #endif // TEST_UFO_METOFFICEBUDDYPAIRFINDER_H_
ufo::test::maxTotalNumBuddies
int maxTotalNumBuddies(const std::vector< ObsPair > &pairs)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:78
ufo::test::CASE
CASE("ufo/MetOfficeBuddyPairFinder/" "Duplicates, constraints on buddy counts, legacy pair collector")
Definition: test/ufo/MetOfficeBuddyPairFinder.h:171
ufo::test::testInvarianceToLongitude
void testInvarianceToLongitude(const eckit::LocalConfiguration &conf)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:191
ufo::test::testDuplicatesAndBuddyCountConstraints
void testDuplicatesAndBuddyCountConstraints(const eckit::LocalConfiguration &conf)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:94
ufo::MetOfficeBuddyPair::rotationBInRad
double rotationBInRad
Definition: MetOfficeBuddyPair.h:25
ufo::test::duplicatePairsPresent
bool duplicatePairsPresent(const std::vector< ObsPair > &pairs)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:37
ufo
Definition: RunCRTM.h:27
ufo::test::MetOfficeBuddyPairFinder::clear
void clear() const override
Definition: test/ufo/MetOfficeBuddyPairFinder.h:354
MetOfficeBuddyPairFinder.h
ufo::Constants::rad2deg
static constexpr double rad2deg
Definition: Constants.h:22
ufo::MetOfficeBuddyPair
Properties of a pair of observations checked against each other during buddy check.
Definition: MetOfficeBuddyPair.h:14
ufo::Constants::mean_earth_rad
static constexpr double mean_earth_rad
Definition: Constants.h:39
ufo::MetOfficeBuddyPair::distanceInKm
double distanceInKm
Definition: MetOfficeBuddyPair.h:23
ufo::MetOfficeBuddyCheckParameters
Options controlling the operation of the MetOfficeBuddyCheck filter.
Definition: MetOfficeBuddyCheckParameters.h:41
ufo::MetOfficeBuddyCheckParameters::searchRadius
oops::Parameter< float > searchRadius
Maximum distance between two observations that may be classified as buddies, in km.
Definition: MetOfficeBuddyCheckParameters.h:49
ufo::MetOfficeBuddyCheckParameters::maxTotalNumBuddies
oops::Parameter< int > maxTotalNumBuddies
Definition: MetOfficeBuddyCheckParameters.h:85
ufo::MetOfficeBuddyPair::rotationAInRad
double rotationAInRad
Definition: MetOfficeBuddyPair.h:24
ufo::test::maxValue
Value maxValue(const std::map< Key, Value > &map)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:67
ufo::MetOfficeBuddyPair::obsIdB
int obsIdB
Definition: MetOfficeBuddyPair.h:22
ufo::test::maxNumBuddiesWithSameStationId
int maxNumBuddiesWithSameStationId(const std::vector< ObsPair > &pairs, const std::vector< int > &stationIds)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:85
ufo::MetOfficeBuddyPair::obsIdA
int obsIdA
Definition: MetOfficeBuddyPair.h:21
ufo::test::contains
bool contains(const std::set< T > &set, const T &element)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:282
ufo::Constants::deg2rad
static constexpr double deg2rad
Definition: Constants.h:21
ufo::test::findBuddyPairs
std::vector< MetOfficeBuddyPair > findBuddyPairs(const MetOfficeBuddyCheckParameters &options, const std::vector< float > &latitudes, const std::vector< float > &longitudes, const std::vector< util::DateTime > &datetimes, const std::vector< float > *pressures, const std::vector< int > &stationIds, const std::vector< size_t > &validObsIds)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:179
ufo::test::reverse
ObsPair reverse(const ObsPair &pair)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:286
ufo::test::MetOfficeBuddyPairFinder
Definition: test/ufo/MetOfficeBuddyPairFinder.h:348
StringUtils.h
ufo::test::testSearchRadius
void testSearchRadius(const eckit::LocalConfiguration &conf)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:290
ufo::test::MetOfficeBuddyPairFinder::register_tests
void register_tests() const override
Definition: test/ufo/MetOfficeBuddyPairFinder.h:352
ufo::test::MetOfficeBuddyPairFinder::testid
std::string testid() const override
Definition: test/ufo/MetOfficeBuddyPairFinder.h:350
ufo::MetOfficeBuddyCheckParameters::maxNumBuddiesWithSameStationId
oops::Parameter< int > maxNumBuddiesWithSameStationId
Definition: MetOfficeBuddyCheckParameters.h:95
ufo::test::ObsPair
std::pair< int, int > ObsPair
Definition: test/ufo/MetOfficeBuddyPairFinder.h:35
MetOfficeBuddyCheckParameters.h
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
conf
Definition: conf.py:1
ufo::test::findEndpoint
void findEndpoint(double startLat, double startLon, double distance, double bearing, double &endLat, double &endLon)
Definition: test/ufo/MetOfficeBuddyPairFinder.h:257