diff --git a/src/ufo/filters/CMakeLists.txt b/src/ufo/filters/CMakeLists.txt index 4382a2107..993bcb480 100644 --- a/src/ufo/filters/CMakeLists.txt +++ b/src/ufo/filters/CMakeLists.txt @@ -78,6 +78,8 @@ set ( filters_files ObsDomainErrCheck.h ObsFilterData.cc ObsFilterData.h + ObsPolygonCheck.cc + 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..e05b0dc74 --- /dev/null +++ b/src/ufo/filters/ObsPolygonCheck.cc @@ -0,0 +1,139 @@ +/* + * 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 +#include +#include +#include + +#include + +#include +#include "ioda/ObsSpace.h" +#include "oops/util/Logger.h" + +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; +} + +// ----------------------------------------------------------------------------- + +ObsPolygonCheck::~ObsPolygonCheck() { + oops::Log::trace() << "ObsPolygonCheck destructor" << std::endl; +} + +// ----------------------------------------------------------------------------- + +void ObsPolygonCheck::applyFilter(const std::vector &apply, + 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 from the parameters a point within the polygon. + 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; + 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 + << ") and vertex latitude count (" << nlat << ")."; + throw eckit::BadValue(what.str(), Here()); + } + 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); + + // 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: boost::geometry does not like your polygon (\"" << reason << "\")"; + throw eckit::BadValue(what.str(), Here()); + } + + // Get the observation locations. + std::vector 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. + // 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" + 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 < nlocs; iloc++) { + if (apply[iloc]) { + 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++) + if (apply[iloc]) + vec[iloc] = notInside[iloc]; + + oops::Log::trace() << "ObsPolygonCheck applyFilter complete (kept " << insideCount + << " of " << applyCount << " locations discarding " + << (applyCount - insideCount) << ")" << 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..153b9d6f5 --- /dev/null +++ b/src/ufo/filters/ObsPolygonCheck.h @@ -0,0 +1,88 @@ +/* + * 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) + + public: + oops::RequiredParameter> vertex_longitudes { + "vertex longitudes", + "Longitudes 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}; +}; + +/// 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 { + 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 diff --git a/test/testinput/unit_tests/filters/CMakeLists.txt b/test/testinput/unit_tests/filters/CMakeLists.txt index 40f071d2c..3259ae2ee 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..7577f8106 --- /dev/null +++ b/test/testinput/unit_tests/filters/polygon_check.yaml @@ -0,0 +1,22 @@ +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 + failedBenchmark: 230