Geographical polygon filter

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

src/ufo/filters/ObsPolygonCheck.cc

@@ -0,0 +1,130 @@
+/*
+ * 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 <boost/geometry.hpp>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace ufo {
+
+// -----------------------------------------------------------------------------
+
+ObsPolygonCheck::ObsPolygonCheck(ioda::ObsSpace & obsdb, const Parameters_ & parameters,
+                                 std::shared_ptr<ioda::ObsDataVector<int> > flags,
+                                 std::shared_ptr<ioda::ObsDataVector<float> > 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 std::vector<bool> &apply,
+                                  const Variables & filtervars,
+                                  std::vector<std::vector<bool>> & 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<double, 2, bg::cs::geographic<bg::degree>>;
+
+  // polygon_t = a closed polygon of those points
+  using polygon_t = bg::model::polygon<point_t>;
+
+  // 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<nlon; i++)
+    poly.outer().push_back(point_t(vertex_longitudes[i], vertex_latitudes[i]));
+
+  // 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)) {
+    std::ostringstream what;
+    what << "ObsPolygonCheck: unable to correct invalid polygon (" << reason << ")";
+    throw ObsPolygonIsInvalid(what.str());
+  }
+
+  // Get the observation locations.
+  std::vector <float> 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.
+  bool useThisSide = bg::within(insidePoint, poly);
+
+  // Find all points that are on the opposite side from the "inside point"
+  std::vector<bool> notInside(obsdb_.nlocs(), true);
+  size_t applyCount = 0;
+  size_t insideCount = 0;
+  for (size_t iloc = 0; iloc < obsdb_.nlocs(); iloc++)
+    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)
+    for (size_t iloc = 0; iloc < obsdb_.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

src/ufo/filters/ObsPolygonCheck.h

@@ -0,0 +1,104 @@
+/*
+ * 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 <exception>
+#include <memory>
+#include <ostream>
+#include <string>
+#include <vector>
+
+#include "oops/util/ObjectCounter.h"
+#include "ufo/filters/FilterBase.h"
+#include "ufo/filters/QCflags.h"
+
+namespace ioda {
+  template <typename DATATYPE> class ObsDataVector;
+  class ObsSpace;
+}
+
+namespace ufo {
+
+class ObsPolygonCheckParameters : public FilterParametersBase {
+  OOPS_CONCRETE_PARAMETERS(ObsPolygonCheckParameters, FilterParametersBase)
+
+ public:
+  oops::RequiredParameter<std::vector<float>> vertex_longitudes {
+    "vertex longitudes",
+    "Longitudes of vertices of the polygon.",
+    this};
+  oops::RequiredParameter<std::vector<float>> vertex_latitudes {
+    "vertex latitudes",
+    "Latitudes of vertices of the polygon.",
+    this};
+  oops::RequiredParameter<float> inside_point_longitude {
+    "inside point longitude",
+    "Longitude of a point inside the polygon (used to determine which side is inside).",
+    this};
+  oops::RequiredParameter<float> 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.
+
+class ObsPolygonCheck : public FilterBase,
+                       private util::ObjectCounter<ObsPolygonCheck> {
+ 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<ioda::ObsDataVector<int> >,
+                 std::shared_ptr<ioda::ObsDataVector<float> >);
+  ~ObsPolygonCheck();
+
+ private:
+  void print(std::ostream &) const override;
+  void applyFilter(const std::vector<bool> &, const Variables &,
+                   std::vector<std::vector<bool>> &) const override;
+  int qcFlag() const override {return QCflags::domain;}
+
+  Parameters_ parameters_;
+};
+
+}  // namespace ufo
+
+#endif  // UFO_FILTERS_OBSPOLYGONCHECK_H_

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<ObsDomainErrCheck>
            domainErrCheckMaker("DomainErr Check");
+  static FilterMaker<ObsPolygonCheck>
+           polygonCheckMaker("Polygon Check");
   static FilterMaker<FinalCheck>
            finalCheckMaker("Final Check");
   static FilterMaker<Gaussian_Thinning>