8 #ifndef TEST_UFO_METOFFICEBUDDYPAIRFINDER_H_
9 #define TEST_UFO_METOFFICEBUDDYPAIRFINDER_H_
18 #include <boost/optional.hpp>
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"
38 bool duplicatesFound =
false;
39 std::set<ObsPair> pairSet;
42 for (
const ObsPair &pair : pairs) {
43 bool inserted = pairSet.insert(pair).second;
45 duplicatesFound =
true;
46 oops::Log::trace() <<
"Duplicate pair found: (" << pair.first <<
", " << pair.second <<
")\n";
51 for (
const ObsPair &pair : pairs) {
52 ObsPair reversePair(pair.second, pair.first);
53 bool inserted = pairSet.insert(reversePair).second;
55 duplicatesFound =
true;
56 oops::Log::trace() <<
"Pair (" << pair.first <<
", " << pair.second
57 <<
") present also in reverse order\n";
61 oops::Log::trace().flush();
63 return duplicatesFound;
66 template <
typename Key,
typename Value>
67 Value
maxValue(
const std::map<Key, Value> &map)
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;
79 std::map<int, int> numBuddiesByObsId;
80 for (
const ObsPair & pair : pairs)
81 ++numBuddiesByObsId[pair.first];
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);
95 util::DateTime bgn(
conf.getString(
"window begin"));
96 util::DateTime end(
conf.getString(
"window end"));
98 const eckit::LocalConfiguration obsSpaceConf(
conf,
"obs space");
99 ioda::ObsSpace obsSpace(obsSpaceConf, oops::mpi::world(), bgn, end, oops::mpi::myself());
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);
107 std::vector<float> latitudes(obsSpace.nlocs());
108 obsSpace.get_db(
"MetaData",
"latitude", latitudes);
110 std::vector<float> longitudes(obsSpace.nlocs());
111 obsSpace.get_db(
"MetaData",
"longitude", longitudes);
113 std::vector<util::DateTime> datetimes(obsSpace.nlocs());
114 obsSpace.get_db(
"MetaData",
"datetime", datetimes);
116 std::vector<int> stationIds(obsSpace.recnum().begin(), obsSpace.recnum().end());
118 std::vector<size_t> validObsIds;
119 if (
conf.has(
"valid_obs_ids")) {
120 validObsIds =
conf.getUnsignedVector(
"valid_obs_ids");
122 validObsIds.resize(obsSpace.nlocs());
123 std::iota(validObsIds.begin(), validObsIds.end(), 0);
126 const eckit::LocalConfiguration modernFilterConf(
conf,
"Met Office Buddy Check, modern");
128 modernOptions.deserialize(modernFilterConf);
130 std::vector<ObsPair> modernPairs;
133 airPressures.get_ptr(), stationIds);
134 const std::vector<MetOfficeBuddyPair> buddyPairs = finder.findBuddyPairs(validObsIds);
137 modernPairs.push_back(
ObsPair(pair.obsIdA, pair.obsIdB));
146 const eckit::LocalConfiguration legacyFilterConf(
conf,
"Met Office Buddy Check, legacy");
148 legacyOptions.deserialize(legacyFilterConf);
150 std::vector<ObsPair> legacyPairs;
153 airPressures.get_ptr(), stationIds);
154 const std::vector<MetOfficeBuddyPair> buddyPairs = finder.findBuddyPairs(validObsIds);
157 legacyPairs.push_back(
ObsPair(pair.obsIdA, pair.obsIdB));
168 EXPECT(modernPairs.size() > legacyPairs.size());
171 CASE(
"ufo/MetOfficeBuddyPairFinder/"
172 "Duplicates, constraints on buddy counts, legacy pair collector") {
174 ::test::TestEnvironment::config(),
175 "Duplicates, constraints on buddy counts, "
176 "legacy pair collector"));
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) {
187 pressures, stationIds);
188 return finder.findBuddyPairs(validObsIds);
192 const float searchRadius = 100;
194 std::vector<float> referenceLatitudes =
conf.getFloatVector(
"reference.latitudes");
195 std::vector<float> referenceLongitudes =
conf.getFloatVector(
"reference.longitudes");
197 std::vector<util::DateTime> datetimes(2, util::DateTime());
198 std::vector<int> stationIds{1, 2};
200 std::vector<size_t> validObsIds{0, 1};
202 const eckit::LocalConfiguration filterConf(
conf,
"Met Office Buddy Check");
204 options.deserialize(filterConf);
206 const std::vector<MetOfficeBuddyPair> referenceBuddyPairs =
findBuddyPairs(
207 options, referenceLatitudes, referenceLongitudes,
208 datetimes,
nullptr, stationIds, validObsIds);
209 EXPECT_EQUAL(referenceBuddyPairs.size(), 1);
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);
224 referenceLatitudes[refPair.
obsIdB] - referenceLatitudes[refPair.
obsIdA] &&
226 referenceLongitudes[refPair.
obsIdB] - referenceLongitudes[refPair.
obsIdA]) {
235 referenceLatitudes[refPair.
obsIdA] - referenceLatitudes[refPair.
obsIdB] &&
237 referenceLongitudes[refPair.
obsIdA] - referenceLongitudes[refPair.
obsIdB]) {
245 CASE(
"ufo/TemporalThinning/Invariance to longitude, different zonal bands") {
247 ::test::TestEnvironment::config(),
248 "Invariance to longitude, different zonal bands"));
251 CASE(
"ufo/TemporalThinning/Invariance to longitude, same zonal band") {
253 ::test::TestEnvironment::config(),
254 "Invariance to longitude, same zonal band"));
258 double &endLat,
double &endLon)
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));
273 float &endLat,
float &endLon)
275 double dEndLat, dEndLon;
281 template <
typename T>
282 bool contains(
const std::set<T> &set,
const T &element) {
283 return set.find(element) != set.end();
287 return ObsPair(pair.second, pair.first);
291 const eckit::LocalConfiguration filterConf(
conf,
"Met Office Buddy Check");
293 options.deserialize(filterConf);
297 for (
const eckit::LocalConfiguration &testConf :
conf.getSubConfigurations(
"test")) {
298 float centralLatitude = testConf.getFloat(
"central_latitude");
299 float centralLongitude = testConf.getFloat(
"central_longitude");
301 std::vector<float> latitudes{centralLatitude};
302 std::vector<float> longitudes{centralLongitude};
303 std::vector<ObsPair> pairsExpectedToBePresent, pairsExpectedToBeAbsent;
305 const int numBearings = 8;
306 for (
float i = 0; i < numBearings; ++i) {
307 float latitude, longitude;
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));
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));
322 const std::vector<util::DateTime> datetimes(latitudes.size(), util::DateTime());
323 const std::vector<int> stationIds(latitudes.size(), 0);
325 std::vector<size_t> validObsIds(latitudes.size());
326 std::iota(validObsIds.begin(), validObsIds.end(), 0);
328 const std::vector<MetOfficeBuddyPair> buddyPairs =
findBuddyPairs(
329 options, latitudes, longitudes, datetimes,
nullptr, stationIds, validObsIds);
331 std::set<ObsPair> pairs;
333 pairs.insert(
ObsPair(pair.obsIdA, pair.obsIdB));
336 for (
const ObsPair &pair : pairsExpectedToBePresent)
339 for (
const ObsPair &pair : pairsExpectedToBeAbsent)
344 CASE(
"ufo/TemporalThinning/Search radius") {
345 testSearchRadius(eckit::LocalConfiguration(::test::TestEnvironment::config(),
"Search radius"));
350 std::string
testid()
const override {
return "ufo::test::MetOfficeBuddyPairFinder";}
360 #endif // TEST_UFO_METOFFICEBUDDYPAIRFINDER_H_