From 9024622a4d2a6f198c19f04ea3e961211d40ba6c Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Wed, 14 Jan 2026 13:45:44 -0500 Subject: [PATCH 01/12] filter for obs within a lon-lat polygon --- src/ufo/filters/CMakeLists.txt | 1 + src/ufo/filters/ObsPolygonCheck.cc | 114 ++++++++++++++++++++++++++ src/ufo/filters/ObsPolygonCheck.h | 80 ++++++++++++++++++ src/ufo/instantiateObsFilterFactory.h | 3 + 4 files changed, 198 insertions(+) create mode 100644 src/ufo/filters/ObsPolygonCheck.cc create mode 100644 src/ufo/filters/ObsPolygonCheck.h diff --git a/src/ufo/filters/CMakeLists.txt b/src/ufo/filters/CMakeLists.txt index 4382a2107..ef7cd4303 100644 --- a/src/ufo/filters/CMakeLists.txt +++ b/src/ufo/filters/CMakeLists.txt @@ -78,6 +78,7 @@ set ( filters_files ObsDomainErrCheck.h ObsFilterData.cc ObsFilterData.h + ObsPolygonCheck.h ObsRefractivityGradientCheck.h ObsRefractivityGradientCheck.cc ParameterSubstitution.cc diff --git a/src/ufo/filters/ObsPolygonCheck.cc b/src/ufo/filters/ObsPolygonCheck.cc new file mode 100644 index 000000000..615afcbd4 --- /dev/null +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -0,0 +1,114 @@ +/* + * CC0 - No Copyright (Public Domain) + * + * The person who associated a work with this deed has dedicated the work to the public domain by waiving all of his or her rights to the work worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law. + * + * You can copy, modify, distribute and perform the work, even for commercial purposes, all without asking permission. See Other Information below. + * + * Other Information + * In no way are the patent or trademark rights of any person affected by CC0, nor are the rights that other persons may have in the work or in how the work is used, such as publicity or privacy rights. + * + * The person who associated a work with this deed makes no warranties about the work, and disclaims liability for all uses of the work, to the fullest extent permitted by applicable law. + * + * When using or citing the work, you should not imply endorsement by the author or the affirmer. + */ + +#include "ufo/filters/ObsPolygonCheck.h" + +#include "ioda/ObsSpace.h" +#include "oops/util/Logger.h" +#include +#include +#include +#include +#include +#include + +namespace { + // Read the entire contents of a text file and return it as a string. + std::string readWholeFile(const std::string &filename) { + return (std::ostringstream() << std::ifstream(filename).rdbuf()).str(); + } +} + +namespace ufo { + +// ----------------------------------------------------------------------------- + +ObsPolygonCheck::ObsPolygonCheck(ioda::ObsSpace & obsdb, const Parameters_ & parameters, + std::shared_ptr > flags, + std::shared_ptr > obserr) + : FilterBase(obsdb, parameters, flags, obserr), parameters_(parameters) +{ + oops::Log::trace() << "ObsPolygonCheck constructor" << std::endl; + oops::Log::debug() << "ObsPolygonCheck: config = " << parameters_ << std::endl; +} + +// ----------------------------------------------------------------------------- + +ObsPolygonCheck::~ObsPolygonCheck() { + oops::Log::trace() << "ObsPolygonCheck destructor" << std::endl; +} + +// ----------------------------------------------------------------------------- + +void ObsPolygonCheck::applyFilter(const Variables & filtervars, + std::vector> & flagged) const { + oops::Log::trace() << "ObsPolygonCheck applyFilter start" << std::endl; + + namespace bg = boost::geometry; + + // point_t = a boost::geometry point as a longitude-latitude pair in degrees + using point_t = bg::model::point>; + + // polygon_t = a closed polygon of those points + using polygon_t = bg::model::polygon; + + // Get the observation locations. + ioda::ObsDataVector lons(obsdb_, "longitude", "MetaData"); + ioda::ObsDataVector lats(obsdb_, "latitude", "MetaData"); + + // Read the polygon from a wkt-format file. + polygon_t poly = bg::from_wkt(readWholeFile(filename)); + + // Ask boost::geometry to correct common problems in the polygon definition. + bg::correct(poly); + + // Scan for obvious errors that bg::correct couldn't correct. + if(std::string reason; !bg::is_valid(poly, reason)) { + oops::Log::error() << "ObsPolygonCheck: unable to correct invalid polygon (" + << reason << "): " << bg::wkt(poly) << std::endl; + return; + } + + // Figure out which side is the inside by checking a point that is known to be inside. + point_t insidePoint = point_t(parmeters_.insideLon, parameters_.insideLat); + bool useThisSide = bg::within(insidePoint, poly); + + // Find all points that are on the opposite side from the "insidePoint" + vector notInside(obsdb_.nlocs(), true); + for (size_t iloc = 0; iloc < obsdb_.nlocs(); iloc++) + try { + bool inside = useThisSide == bg::within(point_t(lons[iloc], lats[iloc]), poly); + notInside[iloc] = not inside; + } catch(const bg::exception &ex) { + // We only catch boost geometry exceptions here; anything else is passed through. + oops::Log::error() << "ObsPolygonCheck: boost::geometry error: " << ex.what() << std::endl; + // The default value for all points is "not inside" so any erroring points will be rejected. + } + + for (auto &vec : flagged) + vec.assign(notInside.begin(), notInside.end()); + + oops::Log::trace() << "ObsPolygonCheck applyFilter complete" << std::endl; +} + +// ----------------------------------------------------------------------------- + +void ObsPolygonCheck::print(std::ostream & os) const { + os << "ObsPolygonCheck: config = " << parameters_ << std::endl; +} + +// ----------------------------------------------------------------------------- + +} // namespace ufo diff --git a/src/ufo/filters/ObsPolygonCheck.h b/src/ufo/filters/ObsPolygonCheck.h new file mode 100644 index 000000000..144c4dbb3 --- /dev/null +++ b/src/ufo/filters/ObsPolygonCheck.h @@ -0,0 +1,80 @@ +/* + * CC0 - No Copyright (Public Domain) + * + * The person who associated a work with this deed has dedicated the work to the public domain by waiving all of his or her rights to the work worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law. + * + * You can copy, modify, distribute and perform the work, even for commercial purposes, all without asking permission. See Other Information below. + * + * Other Information + * In no way are the patent or trademark rights of any person affected by CC0, nor are the rights that other persons may have in the work or in how the work is used, such as publicity or privacy rights. + * + * The person who associated a work with this deed makes no warranties about the work, and disclaims liability for all uses of the work, to the fullest extent permitted by applicable law. + * + * When using or citing the work, you should not imply endorsement by the author or the affirmer. + */ + +#ifndef UFO_FILTERS_OBSPOLYGONCHECK_H_ +#define UFO_FILTERS_OBSPOLYGONCHECK_H_ + +#include +#include +#include +#include + +#include "oops/util/ObjectCounter.h" +#include "ufo/filters/FilterBase.h" +#include "ufo/filters/QCflags.h" + +namespace ioda { + template class ObsDataVector; + class ObsSpace; +} + +namespace ufo { + +class ObsPolygonCheckParameters : public FilterParametersBase { + OOPS_CONCRETE_PARAMETERS(ObsPolygonCheckParameters, FilterParametersBase) + oops::RequiredParameter polyFileName + {"polygon wkt file", + "Path to a Well Known Text (WKT) file containing the polygon as longitude, latitude pairs. " + "Observations that lie within this polygon are retained. " + "File contents should look like: \"POLYGON((180 -3, -125 18, 98.4 0, 180 -3))\"" + this}; + oops::RequiredParameter insideLon + {"inside point longitude", + "Longitude of the sample point within the domain, used to find which side of the polygon is inside.", + this}; + oops::RequiredParameter insideLat + {"inside point latitude", + "Latitude of the sample point within the domain, used to find which side of the polygon is inside.", + this}; +}; + +/// PolygonCheck: find obs within a specified polygon. + +class ObsPolygonCheck : public FilterBase, + private util::ObjectCounter { + public: + /// The type of parameters accepted by the constructor of this filter. + /// This typedef is used by the FilterFactory. + typedef ObsPolygonCheckParameters Parameters_; + + static const std::string classname() {return "ufo::ObsPolygonCheck";} + + ObsPolygonCheck(ioda::ObsSpace &, const Parameters_ &, + std::shared_ptr >, + std::shared_ptr >); + ~ObsPolygonCheck(); + + private: + void print(std::ostream &) const override; + void applyFilter(const std::vector &, const Variables &, + std::vector> &) const override; + int qcFlag() const override {return QCflags::domain;} + + Parameters_ parameters_; +}; + +} // namespace ufo + +#endif // UFO_FILTERS_OBSPOLYGONCHECK_H_ diff --git a/src/ufo/instantiateObsFilterFactory.h b/src/ufo/instantiateObsFilterFactory.h index c9c3798e7..15816398b 100644 --- a/src/ufo/instantiateObsFilterFactory.h +++ b/src/ufo/instantiateObsFilterFactory.h @@ -36,6 +36,7 @@ #include "ufo/filters/ObsDiagnosticsWriter.h" #include "ufo/filters/ObsDomainCheck.h" #include "ufo/filters/ObsDomainErrCheck.h" +#include "ufo/filters/ObsPolygonCheck.h" #include "ufo/filters/ObsRefractivityGradientCheck.h" #include "ufo/filters/ParameterSubstitution.h" #include "ufo/filters/PerformAction.h" @@ -110,6 +111,8 @@ void instantiateObsFilterFactory() { domainCheckMaker("Domain Check"); static FilterMaker domainErrCheckMaker("DomainErr Check"); + static FilterMaker + polygonCheckMaker("Polygon Check"); static FilterMaker finalCheckMaker("Final Check"); static FilterMaker From bab8034e9c6e06215bc9f957f93d4bd5d369e58a Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Tue, 27 Jan 2026 18:44:19 -0500 Subject: [PATCH 02/12] add ObsPolygonCheck to CMakeLists --- src/ufo/filters/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ufo/filters/CMakeLists.txt b/src/ufo/filters/CMakeLists.txt index ef7cd4303..993bcb480 100644 --- a/src/ufo/filters/CMakeLists.txt +++ b/src/ufo/filters/CMakeLists.txt @@ -78,6 +78,7 @@ set ( filters_files ObsDomainErrCheck.h ObsFilterData.cc ObsFilterData.h + ObsPolygonCheck.cc ObsPolygonCheck.h ObsRefractivityGradientCheck.h ObsRefractivityGradientCheck.cc From fe2c5c89faa122ef688e1a678152b485a108e9bb Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Tue, 27 Jan 2026 18:44:30 -0500 Subject: [PATCH 03/12] Move polygon into yaml --- src/ufo/filters/ObsPolygonCheck.cc | 78 ++++++++++++++++++------------ src/ufo/filters/ObsPolygonCheck.h | 52 ++++++++++++++------ 2 files changed, 85 insertions(+), 45 deletions(-) diff --git a/src/ufo/filters/ObsPolygonCheck.cc b/src/ufo/filters/ObsPolygonCheck.cc index 615afcbd4..6c052bfae 100644 --- a/src/ufo/filters/ObsPolygonCheck.cc +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -18,19 +18,11 @@ #include "ioda/ObsSpace.h" #include "oops/util/Logger.h" #include -#include #include #include #include #include -namespace { - // Read the entire contents of a text file and return it as a string. - std::string readWholeFile(const std::string &filename) { - return (std::ostringstream() << std::ifstream(filename).rdbuf()).str(); - } -} - namespace ufo { // ----------------------------------------------------------------------------- @@ -52,7 +44,8 @@ ObsPolygonCheck::~ObsPolygonCheck() { // ----------------------------------------------------------------------------- -void ObsPolygonCheck::applyFilter(const Variables & filtervars, +void ObsPolygonCheck::applyFilter(const std::vector &apply, + const Variables & filtervars, std::vector> & flagged) const { oops::Log::trace() << "ObsPolygonCheck applyFilter start" << std::endl; @@ -64,43 +57,66 @@ void ObsPolygonCheck::applyFilter(const Variables & filtervars, // polygon_t = a closed polygon of those points using polygon_t = bg::model::polygon; - // Get the observation locations. - ioda::ObsDataVector lons(obsdb_, "longitude", "MetaData"); - ioda::ObsDataVector lats(obsdb_, "latitude", "MetaData"); - - // Read the polygon from a wkt-format file. - polygon_t poly = bg::from_wkt(readWholeFile(filename)); + // Get from the parameters a point within the polygon. + point_t insidePoint(parameters_.inside_point_longitude.value(), parameters_.inside_point_latitude.value()); + + // Assemble a polygon from the list of longitudes and latitudes. + polygon_t poly; + auto &vertex_latitudes = parameters_.vertex_latitudes.value(); + auto &vertex_longitudes = parameters_.vertex_longitudes.value(); + size_t nlon = vertex_longitudes.size(); + size_t nlat = vertex_latitudes.size(); + if(nlon != nlat) { + std::ostringstream what; + what << "Mismatch between vertex longitude count (" << nlon + << ") and vertex latitude count (" << nlat << ")."; + throw ObsPolygonLatLonSizeMismatch(what.str()); + } + for(int i=0; i lats, lons; + obsdb_.get_db("MetaData", "latitude", lats); + obsdb_.get_db("MetaData", "longitude", lons); + // Figure out which side is the inside by checking a point that is known to be inside. - point_t insidePoint = point_t(parmeters_.insideLon, parameters_.insideLat); bool useThisSide = bg::within(insidePoint, poly); - // Find all points that are on the opposite side from the "insidePoint" - vector notInside(obsdb_.nlocs(), true); + // Find all points that are on the opposite side from the "inside point" + std::vector notInside(obsdb_.nlocs(), true); + size_t applyCount = 0; + size_t insideCount = 0; for (size_t iloc = 0; iloc < obsdb_.nlocs(); iloc++) - try { - bool inside = useThisSide == bg::within(point_t(lons[iloc], lats[iloc]), poly); - notInside[iloc] = not inside; - } catch(const bg::exception &ex) { - // We only catch boost geometry exceptions here; anything else is passed through. - oops::Log::error() << "ObsPolygonCheck: boost::geometry error: " << ex.what() << std::endl; - // The default value for all points is "not inside" so any erroring points will be rejected. - } + if(apply[iloc]) + try { + applyCount++; + bool inside = useThisSide == bg::within(point_t(lons[iloc], lats[iloc]), poly); + notInside[iloc] = not inside; + if(inside) + insideCount++; + } catch(const bg::exception &ex) { + // We only catch boost geometry exceptions here; anything else is passed through. + oops::Log::error() << "ObsPolygonCheck: boost::geometry error: " << ex.what() << std::endl; + // The default value for all points is "not inside" so any erroring points will be rejected. + } for (auto &vec : flagged) - vec.assign(notInside.begin(), notInside.end()); + for (size_t iloc = 0; iloc < obsdb_.nlocs(); iloc++) + if(apply[iloc]) + vec[iloc] = notInside[iloc]; - oops::Log::trace() << "ObsPolygonCheck applyFilter complete" << std::endl; + oops::Log::trace() << "ObsPolygonCheck applyFilter complete (kept " << insideCount << " of " << applyCount << " locations discarding " << (applyCount - insideCount) << ")" << std::endl; } // ----------------------------------------------------------------------------- diff --git a/src/ufo/filters/ObsPolygonCheck.h b/src/ufo/filters/ObsPolygonCheck.h index 144c4dbb3..aa74e36b8 100644 --- a/src/ufo/filters/ObsPolygonCheck.h +++ b/src/ufo/filters/ObsPolygonCheck.h @@ -16,6 +16,7 @@ #ifndef UFO_FILTERS_OBSPOLYGONCHECK_H_ #define UFO_FILTERS_OBSPOLYGONCHECK_H_ +#include #include #include #include @@ -34,20 +35,43 @@ namespace ufo { class ObsPolygonCheckParameters : public FilterParametersBase { OOPS_CONCRETE_PARAMETERS(ObsPolygonCheckParameters, FilterParametersBase) - oops::RequiredParameter polyFileName - {"polygon wkt file", - "Path to a Well Known Text (WKT) file containing the polygon as longitude, latitude pairs. " - "Observations that lie within this polygon are retained. " - "File contents should look like: \"POLYGON((180 -3, -125 18, 98.4 0, 180 -3))\"" - this}; - oops::RequiredParameter insideLon - {"inside point longitude", - "Longitude of the sample point within the domain, used to find which side of the polygon is inside.", - this}; - oops::RequiredParameter insideLat - {"inside point latitude", - "Latitude of the sample point within the domain, used to find which side of the polygon is inside.", - this}; + + public: + oops::RequiredParameter> vertex_longitudes { + "vertex longitudes", + "Latitudes of vertices of the polygon.", + this}; + oops::RequiredParameter> vertex_latitudes { + "vertex latitudes", + "Latitudes of vertices of the polygon.", + this}; + oops::RequiredParameter inside_point_longitude { + "inside point longitude", + "Longitude of a point inside the polygon (used to determine which side is inside).", + this}; + oops::RequiredParameter inside_point_latitude { + "inside point latitude", + "Latitude of a point inside the polygon (used to determine which side is inside).", + this}; +}; + +/// ObsPolygonLatLonSizeMismatch: thrown when the parameters vertex_longitudes +/// and vertex_latitudes have different lengths. + +class ObsPolygonLatLonSizeMismatch: public std::invalid_argument { +public: + ObsPolygonLatLonSizeMismatch(const std::string &message): + std::invalid_argument(message) + {} +}; + +/// ObsPolygonIsInvalid: thrown when boost::geometry::is_valid doesn't like a polygon. + +class ObsPolygonIsInvalid: public std::invalid_argument { +public: + ObsPolygonIsInvalid(const std::string &message): + std::invalid_argument(message) + {} }; /// PolygonCheck: find obs within a specified polygon. From a694370e78a899cb709fb1ec8512676a3c324980 Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Wed, 28 Jan 2026 12:24:38 -0500 Subject: [PATCH 04/12] Latitude -> Longitude in a help message --- src/ufo/filters/ObsPolygonCheck.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ufo/filters/ObsPolygonCheck.h b/src/ufo/filters/ObsPolygonCheck.h index aa74e36b8..4cdf7fe7c 100644 --- a/src/ufo/filters/ObsPolygonCheck.h +++ b/src/ufo/filters/ObsPolygonCheck.h @@ -39,7 +39,7 @@ class ObsPolygonCheckParameters : public FilterParametersBase { public: oops::RequiredParameter> vertex_longitudes { "vertex longitudes", - "Latitudes of vertices of the polygon.", + "Longitudes of vertices of the polygon.", this}; oops::RequiredParameter> vertex_latitudes { "vertex latitudes", From fbf2b38f77f3ad634468f0531b1d97c0a587df0f Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Wed, 28 Jan 2026 13:31:05 -0500 Subject: [PATCH 05/12] address reviewer comments --- src/ufo/filters/ObsPolygonCheck.cc | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/ufo/filters/ObsPolygonCheck.cc b/src/ufo/filters/ObsPolygonCheck.cc index 6c052bfae..9571b1c1e 100644 --- a/src/ufo/filters/ObsPolygonCheck.cc +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -58,14 +58,14 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, using polygon_t = bg::model::polygon; // Get from the parameters a point within the polygon. - point_t insidePoint(parameters_.inside_point_longitude.value(), parameters_.inside_point_latitude.value()); + const point_t insidePoint(parameters_.inside_point_longitude.value(), parameters_.inside_point_latitude.value()); // Assemble a polygon from the list of longitudes and latitudes. polygon_t poly; - auto &vertex_latitudes = parameters_.vertex_latitudes.value(); - auto &vertex_longitudes = parameters_.vertex_longitudes.value(); - size_t nlon = vertex_longitudes.size(); - size_t nlat = vertex_latitudes.size(); + const auto &vertex_latitudes = parameters_.vertex_latitudes.value(); + const auto &vertex_longitudes = parameters_.vertex_longitudes.value(); + const size_t nlon = vertex_longitudes.size(); + const size_t nlat = vertex_latitudes.size(); if(nlon != nlat) { std::ostringstream what; what << "Mismatch between vertex longitude count (" << nlon @@ -91,13 +91,14 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, obsdb_.get_db("MetaData", "longitude", lons); // Figure out which side is the inside by checking a point that is known to be inside. - bool useThisSide = bg::within(insidePoint, poly); + const bool useThisSide = bg::within(insidePoint, poly); // Find all points that are on the opposite side from the "inside point" - std::vector notInside(obsdb_.nlocs(), true); + const size_t nlocs = obsdb_.nlocs(); + std::vector notInside(nlocs, true); size_t applyCount = 0; size_t insideCount = 0; - for (size_t iloc = 0; iloc < obsdb_.nlocs(); iloc++) + for (size_t iloc = 0; iloc < nlocs; iloc++) if(apply[iloc]) try { applyCount++; @@ -112,7 +113,7 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, } for (auto &vec : flagged) - for (size_t iloc = 0; iloc < obsdb_.nlocs(); iloc++) + for (size_t iloc = 0; iloc < nlocs; iloc++) if(apply[iloc]) vec[iloc] = notInside[iloc]; From c01a48d433cc6fb1390bb8a01eb5a68feb6d04cf Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Thu, 29 Jan 2026 16:54:51 -0500 Subject: [PATCH 06/12] parameters_ are too long to print even as debug due to 12000+ elements --- src/ufo/filters/ObsPolygonCheck.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ufo/filters/ObsPolygonCheck.cc b/src/ufo/filters/ObsPolygonCheck.cc index 9571b1c1e..056d9f3e6 100644 --- a/src/ufo/filters/ObsPolygonCheck.cc +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -33,7 +33,6 @@ ObsPolygonCheck::ObsPolygonCheck(ioda::ObsSpace & obsdb, const Parameters_ & par : FilterBase(obsdb, parameters, flags, obserr), parameters_(parameters) { oops::Log::trace() << "ObsPolygonCheck constructor" << std::endl; - oops::Log::debug() << "ObsPolygonCheck: config = " << parameters_ << std::endl; } // ----------------------------------------------------------------------------- From dba3768970095de01329acd1429620de5ad0e638 Mon Sep 17 00:00:00 2001 From: "Samuel.Trahan" Date: Thu, 29 Jan 2026 21:58:12 +0000 Subject: [PATCH 07/12] meet coding standards --- src/ufo/filters/ObsPolygonCheck.cc | 28 ++++++++++++++++------------ src/ufo/filters/ObsPolygonCheck.h | 8 ++++---- 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/src/ufo/filters/ObsPolygonCheck.cc b/src/ufo/filters/ObsPolygonCheck.cc index 056d9f3e6..d2be390a0 100644 --- a/src/ufo/filters/ObsPolygonCheck.cc +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -15,13 +15,13 @@ #include "ufo/filters/ObsPolygonCheck.h" -#include "ioda/ObsSpace.h" -#include "oops/util/Logger.h" -#include #include #include #include #include +#include +#include "ioda/ObsSpace.h" +#include "oops/util/Logger.h" namespace ufo { @@ -57,7 +57,8 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, using polygon_t = bg::model::polygon; // Get from the parameters a point within the polygon. - const point_t insidePoint(parameters_.inside_point_longitude.value(), parameters_.inside_point_latitude.value()); + const point_t insidePoint(parameters_.inside_point_longitude.value(), + parameters_.inside_point_latitude.value()); // Assemble a polygon from the list of longitudes and latitudes. polygon_t poly; @@ -65,20 +66,20 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, const auto &vertex_longitudes = parameters_.vertex_longitudes.value(); const size_t nlon = vertex_longitudes.size(); const size_t nlat = vertex_latitudes.size(); - if(nlon != nlat) { + if (nlon != nlat) { std::ostringstream what; what << "Mismatch between vertex longitude count (" << nlon << ") and vertex latitude count (" << nlat << ")."; throw ObsPolygonLatLonSizeMismatch(what.str()); } - for(int i=0; i &apply, size_t applyCount = 0; size_t insideCount = 0; for (size_t iloc = 0; iloc < nlocs; iloc++) - if(apply[iloc]) + if (apply[iloc]) { try { applyCount++; bool inside = useThisSide == bg::within(point_t(lons[iloc], lats[iloc]), poly); - notInside[iloc] = not inside; - if(inside) + notInside[iloc] = !inside; + if (inside) insideCount++; } catch(const bg::exception &ex) { // We only catch boost geometry exceptions here; anything else is passed through. oops::Log::error() << "ObsPolygonCheck: boost::geometry error: " << ex.what() << std::endl; // The default value for all points is "not inside" so any erroring points will be rejected. } + } for (auto &vec : flagged) for (size_t iloc = 0; iloc < nlocs; iloc++) - if(apply[iloc]) + if (apply[iloc]) vec[iloc] = notInside[iloc]; - oops::Log::trace() << "ObsPolygonCheck applyFilter complete (kept " << insideCount << " of " << applyCount << " locations discarding " << (applyCount - insideCount) << ")" << std::endl; + oops::Log::trace() << "ObsPolygonCheck applyFilter complete (kept " << insideCount + << " of " << applyCount << " locations discarding " + << (applyCount - insideCount) << ")" << std::endl; } // ----------------------------------------------------------------------------- diff --git a/src/ufo/filters/ObsPolygonCheck.h b/src/ufo/filters/ObsPolygonCheck.h index 4cdf7fe7c..7bb9c30a4 100644 --- a/src/ufo/filters/ObsPolygonCheck.h +++ b/src/ufo/filters/ObsPolygonCheck.h @@ -59,8 +59,8 @@ class ObsPolygonCheckParameters : public FilterParametersBase { /// and vertex_latitudes have different lengths. class ObsPolygonLatLonSizeMismatch: public std::invalid_argument { -public: - ObsPolygonLatLonSizeMismatch(const std::string &message): + public: + explicit ObsPolygonLatLonSizeMismatch(const std::string &message): std::invalid_argument(message) {} }; @@ -68,8 +68,8 @@ class ObsPolygonLatLonSizeMismatch: public std::invalid_argument { /// ObsPolygonIsInvalid: thrown when boost::geometry::is_valid doesn't like a polygon. class ObsPolygonIsInvalid: public std::invalid_argument { -public: - ObsPolygonIsInvalid(const std::string &message): + public: + explicit ObsPolygonIsInvalid(const std::string &message): std::invalid_argument(message) {} }; From d973167ee333c676255ba195654e01afc7dc35e9 Mon Sep 17 00:00:00 2001 From: "Samuel.Trahan" Date: Fri, 30 Jan 2026 22:01:33 +0000 Subject: [PATCH 08/12] ctest for Polygon Check --- .../unit_tests/filters/CMakeLists.txt | 10 +++++++++ .../unit_tests/filters/polygon_check.yaml | 21 +++++++++++++++++++ 2 files changed, 31 insertions(+) create mode 100644 test/testinput/unit_tests/filters/polygon_check.yaml diff --git a/test/testinput/unit_tests/filters/CMakeLists.txt b/test/testinput/unit_tests/filters/CMakeLists.txt index f1da4095c..3f4561f87 100644 --- a/test/testinput/unit_tests/filters/CMakeLists.txt +++ b/test/testinput/unit_tests/filters/CMakeLists.txt @@ -951,6 +951,16 @@ ufo_add_test( NAME test_parameter_substitution LIBS ufo LABELS filters WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../../../ ) +# +ufo_add_test( NAME test_ufo_polygon_check + TIER 1 + ECBUILD + COMMAND ${CMAKE_BINARY_DIR}/bin/test_ObsFilters.x + ARGS "${CMAKE_CURRENT_SOURCE_DIR}/polygon_check.yaml" + MPI 1 + LIBS ufo + LABELS filters + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../../../ ) # Filters that require CRTM in ObsOperator diff --git a/test/testinput/unit_tests/filters/polygon_check.yaml b/test/testinput/unit_tests/filters/polygon_check.yaml new file mode 100644 index 000000000..efc33b841 --- /dev/null +++ b/test/testinput/unit_tests/filters/polygon_check.yaml @@ -0,0 +1,21 @@ +time window: + begin: 2018-04-14T20:30:00Z + end: 2018-04-15T03:30:00Z + +observations: +- obs space: + name: Radiosonde + obsdatain: + engine: + type: H5File + obsfile: Data/ufo/testinput_tier_1/sondes_obs_2018041500_m.nc4 + simulated variables: [airTemperature, windEastward, windNorthward] + obs filters: + - filter: Polygon Check + inside point longitude: 51 + inside point latitude: 59 + vertex longitudes: [ 50, 50, 200, 200, 125, 50 ] + vertex latitudes: [ 20, 60, 60, 20, 50, 20 ] + action: + name: reduce obs space + passedBenchmark: 30 From 13e9a208e81c8d821bce02e6807e38af285b4fb8 Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Mon, 2 Feb 2026 16:42:27 -0500 Subject: [PATCH 09/12] resolve reviewer comments --- src/ufo/filters/ObsPolygonCheck.cc | 39 ++++++++++++++++-------------- src/ufo/filters/ObsPolygonCheck.h | 21 ++++++++++------ 2 files changed, 35 insertions(+), 25 deletions(-) diff --git a/src/ufo/filters/ObsPolygonCheck.cc b/src/ufo/filters/ObsPolygonCheck.cc index d2be390a0..7ca130d75 100644 --- a/src/ufo/filters/ObsPolygonCheck.cc +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -19,6 +19,7 @@ #include #include #include + #include #include "ioda/ObsSpace.h" #include "oops/util/Logger.h" @@ -70,19 +71,23 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, std::ostringstream what; what << "Mismatch between vertex longitude count (" << nlon << ") and vertex latitude count (" << nlat << ")."; - throw ObsPolygonLatLonSizeMismatch(what.str()); + throw ObsPolygonLatLonSizeMismatch(what.str(), Here()); } - for (int i = 0; i < nlon; i++) - poly.outer().push_back(point_t(vertex_longitudes[i], vertex_latitudes[i])); + poly.outer().reserve(nlon); + for (size_t i = 0; i < nlon; i++) + poly.outer().emplace_back(vertex_longitudes[i], vertex_latitudes[i]); // Ask boost::geometry to correct common problems in the polygon definition. + // As of boost 1.91, the only practical effect is to close open rings. + // It will also correct the vertex ordering, but vertex ordering on a + // sphere is ignored. bg::correct(poly); - // Scan for obvious errors that bg::correct couldn't correct. + // Check for polygons that the boost::geometry library doesn't know how to handle. if (std::string reason; !bg::is_valid(poly, reason)) { std::ostringstream what; - what << "ObsPolygonCheck: unable to correct invalid polygon (" << reason << ")"; - throw ObsPolygonIsInvalid(what.str()); + what << "ObsPolygonCheck: boost::geometry does not like your polygon (\"" << reason << "\")"; + throw ObsPolygonIsInvalid(what.str(), Here()); } // Get the observation locations. @@ -91,6 +96,9 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, obsdb_.get_db("MetaData", "longitude", lons); // Figure out which side is the inside by checking a point that is known to be inside. + // As of version 1.91, boost::geometry ignores the clockwise vs. counterclockwise when + // deciding which side of a polygon is the inside in spherical geometry. Hence, + // we need this "inside point" to determine which side is the inside. const bool useThisSide = bg::within(insidePoint, poly); // Find all points that are on the opposite side from the "inside point" @@ -98,20 +106,15 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, std::vector notInside(nlocs, true); size_t applyCount = 0; size_t insideCount = 0; - for (size_t iloc = 0; iloc < nlocs; iloc++) + for (size_t iloc = 0; iloc < nlocs; iloc++) { if (apply[iloc]) { - try { - applyCount++; - bool inside = useThisSide == bg::within(point_t(lons[iloc], lats[iloc]), poly); - notInside[iloc] = !inside; - if (inside) - insideCount++; - } catch(const bg::exception &ex) { - // We only catch boost geometry exceptions here; anything else is passed through. - oops::Log::error() << "ObsPolygonCheck: boost::geometry error: " << ex.what() << std::endl; - // The default value for all points is "not inside" so any erroring points will be rejected. - } + applyCount++; + const bool inside = useThisSide == bg::within(point_t(lons[iloc], lats[iloc]), poly); + notInside[iloc] = !inside; + if (inside) + insideCount++; } + } for (auto &vec : flagged) for (size_t iloc = 0; iloc < nlocs; iloc++) diff --git a/src/ufo/filters/ObsPolygonCheck.h b/src/ufo/filters/ObsPolygonCheck.h index 7bb9c30a4..b78409e2f 100644 --- a/src/ufo/filters/ObsPolygonCheck.h +++ b/src/ufo/filters/ObsPolygonCheck.h @@ -22,6 +22,7 @@ #include #include +#include #include "oops/util/ObjectCounter.h" #include "ufo/filters/FilterBase.h" #include "ufo/filters/QCflags.h" @@ -58,23 +59,29 @@ class ObsPolygonCheckParameters : public FilterParametersBase { /// ObsPolygonLatLonSizeMismatch: thrown when the parameters vertex_longitudes /// and vertex_latitudes have different lengths. -class ObsPolygonLatLonSizeMismatch: public std::invalid_argument { +class ObsPolygonLatLonSizeMismatch: public eckit::BadValue { public: - explicit ObsPolygonLatLonSizeMismatch(const std::string &message): - std::invalid_argument(message) + explicit ObsPolygonLatLonSizeMismatch(const std::string &message, + const eckit::CodeLocation &location = {}): + eckit::BadValue(message, location) {} }; /// ObsPolygonIsInvalid: thrown when boost::geometry::is_valid doesn't like a polygon. -class ObsPolygonIsInvalid: public std::invalid_argument { +class ObsPolygonIsInvalid: public eckit::BadValue { public: - explicit ObsPolygonIsInvalid(const std::string &message): - std::invalid_argument(message) + explicit ObsPolygonIsInvalid(const std::string &message, + const eckit::CodeLocation &location = {}): + eckit::BadValue(message, location) {} }; -/// PolygonCheck: find obs within a specified polygon. +/// PolygonCheck: flags all obs that aren't inside a specified +/// polygon; uses the boost::geometry library. The boost::geometry +/// library doesn't know which side is the inside of the polygon when +/// using spherical geometry. Hence, the caller must provide an +/// "inside point." Points on the other side of the polygon are flagged. class ObsPolygonCheck : public FilterBase, private util::ObjectCounter { From 980e88bff314e95e5bb8f9b8ea57ab670b663cac Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Thu, 5 Feb 2026 11:53:27 -0500 Subject: [PATCH 10/12] remove custom exception classes --- src/ufo/filters/ObsPolygonCheck.cc | 5 +++-- src/ufo/filters/ObsPolygonCheck.h | 23 ----------------------- 2 files changed, 3 insertions(+), 25 deletions(-) diff --git a/src/ufo/filters/ObsPolygonCheck.cc b/src/ufo/filters/ObsPolygonCheck.cc index 7ca130d75..db43a8bf6 100644 --- a/src/ufo/filters/ObsPolygonCheck.cc +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -21,6 +21,7 @@ #include #include +#include #include "ioda/ObsSpace.h" #include "oops/util/Logger.h" @@ -71,7 +72,7 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, std::ostringstream what; what << "Mismatch between vertex longitude count (" << nlon << ") and vertex latitude count (" << nlat << ")."; - throw ObsPolygonLatLonSizeMismatch(what.str(), Here()); + throw eckit::BadValue(what.str(), Here()); } poly.outer().reserve(nlon); for (size_t i = 0; i < nlon; i++) @@ -87,7 +88,7 @@ void ObsPolygonCheck::applyFilter(const std::vector &apply, if (std::string reason; !bg::is_valid(poly, reason)) { std::ostringstream what; what << "ObsPolygonCheck: boost::geometry does not like your polygon (\"" << reason << "\")"; - throw ObsPolygonIsInvalid(what.str(), Here()); + throw eckit::BadValue(what.str(), Here()); } // Get the observation locations. diff --git a/src/ufo/filters/ObsPolygonCheck.h b/src/ufo/filters/ObsPolygonCheck.h index b78409e2f..153b9d6f5 100644 --- a/src/ufo/filters/ObsPolygonCheck.h +++ b/src/ufo/filters/ObsPolygonCheck.h @@ -16,13 +16,11 @@ #ifndef UFO_FILTERS_OBSPOLYGONCHECK_H_ #define UFO_FILTERS_OBSPOLYGONCHECK_H_ -#include #include #include #include #include -#include #include "oops/util/ObjectCounter.h" #include "ufo/filters/FilterBase.h" #include "ufo/filters/QCflags.h" @@ -56,27 +54,6 @@ class ObsPolygonCheckParameters : public FilterParametersBase { this}; }; -/// ObsPolygonLatLonSizeMismatch: thrown when the parameters vertex_longitudes -/// and vertex_latitudes have different lengths. - -class ObsPolygonLatLonSizeMismatch: public eckit::BadValue { - public: - explicit ObsPolygonLatLonSizeMismatch(const std::string &message, - const eckit::CodeLocation &location = {}): - eckit::BadValue(message, location) - {} -}; - -/// ObsPolygonIsInvalid: thrown when boost::geometry::is_valid doesn't like a polygon. - -class ObsPolygonIsInvalid: public eckit::BadValue { - public: - explicit ObsPolygonIsInvalid(const std::string &message, - const eckit::CodeLocation &location = {}): - eckit::BadValue(message, location) - {} -}; - /// PolygonCheck: flags all obs that aren't inside a specified /// polygon; uses the boost::geometry library. The boost::geometry /// library doesn't know which side is the inside of the polygon when From 83566a183ce20c0e7e829c80ae3c92dbfdbf5eb9 Mon Sep 17 00:00:00 2001 From: "Samuel.Trahan" Date: Thu, 5 Feb 2026 19:08:26 +0000 Subject: [PATCH 11/12] add failedBenchmark in Polygon Check unit test --- test/testinput/unit_tests/filters/polygon_check.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/testinput/unit_tests/filters/polygon_check.yaml b/test/testinput/unit_tests/filters/polygon_check.yaml index efc33b841..7577f8106 100644 --- a/test/testinput/unit_tests/filters/polygon_check.yaml +++ b/test/testinput/unit_tests/filters/polygon_check.yaml @@ -19,3 +19,4 @@ observations: action: name: reduce obs space passedBenchmark: 30 + failedBenchmark: 230 From 1f562cff3d566870d21d583daedd3070908a29b3 Mon Sep 17 00:00:00 2001 From: "Samuel.Trahan" Date: Fri, 6 Feb 2026 16:10:30 +0000 Subject: [PATCH 12/12] move header for coding standards --- src/ufo/filters/ObsPolygonCheck.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ufo/filters/ObsPolygonCheck.cc b/src/ufo/filters/ObsPolygonCheck.cc index db43a8bf6..e05b0dc74 100644 --- a/src/ufo/filters/ObsPolygonCheck.cc +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -20,8 +20,9 @@ #include #include -#include #include + +#include #include "ioda/ObsSpace.h" #include "oops/util/Logger.h"