Determine if a geographic area contains a point

scwx-qt/scwx-qt.cmake

@@ -14,6 +14,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
 find_package(Boost)
 find_package(Fontconfig)
 find_package(geographiclib)
+find_package(geos)
 find_package(GLEW)
 find_package(glm)
 find_package(Python COMPONENTS Interpreter)
@@ -533,6 +534,8 @@ target_link_libraries(scwx-qt PUBLIC Qt${QT_VERSION_MAJOR}::Widgets
                                      $<$<CXX_COMPILER_ID:MSVC>:opengl32>
                                      Fontconfig::Fontconfig
                                      GeographicLib::GeographicLib
+                                     GEOS::geos
+                                     GEOS::geos_cxx_flags
                                      GLEW::GLEW
                                      glm::glm
                                      imgui

scwx-qt/source/scwx/qt/util/geographic_lib.cpp

@@ -1,4 +1,9 @@
 #include <scwx/qt/util/geographic_lib.hpp>
+#include <scwx/util/logger.hpp>
+
+#include <GeographicLib/Gnomonic.hpp>
+#include <geos/algorithm/PointLocation.h>
+#include <geos/geom/CoordinateSequence.h>
 
 namespace scwx
 {
@@ -9,6 +14,9 @@ namespace util
 namespace GeographicLib
 {
 
+static const std::string logPrefix_ = "scwx::qt::util::geographic_lib";
+static const auto        logger_    = scwx::util::Logger::Create(logPrefix_);
+
 const ::GeographicLib::Geodesic& DefaultGeodesic()
 {
    static const ::GeographicLib::Geodesic geodesic_ {
@@ -18,6 +26,60 @@ const ::GeographicLib::Geodesic& DefaultGeodesic()
    return geodesic_;
 }
 
+bool AreaContainsPoint(const std::vector<common::Coordinate>& area,
+                       const common::Coordinate&              point)
+{
+   // Cannot have an area with just two points
+   if (area.size() <= 2 || area.size() == 3 && area.front() == area.back())
+   {
+      return false;
+   }
+
+   ::GeographicLib::Gnomonic      gnomonic {};
+   geos::geom::CoordinateSequence sequence {};
+   double                         x;
+   double                         y;
+   bool                           areaContainsPoint = false;
+
+   // Using a gnomonic projection with the test point as the center
+   // latitude/longitude, the projected test point will be at (0, 0)
+   geos::geom::CoordinateXY zero {};
+
+   // Create the area coordinate sequence using a gnomonic projection
+   for (auto& areaCoordinate : area)
+   {
+      gnomonic.Forward(point.latitude_,
+                       point.longitude_,
+                       areaCoordinate.latitude_,
+                       areaCoordinate.longitude_,
+                       x,
+                       y);
+      sequence.add(x, y);
+   }
+
+   // If the sequence is not a ring, add the first point again for closure
+   if (!sequence.isRing())
+   {
+      sequence.add(sequence.front(), false);
+   }
+
+   // The sequence should be a ring at this point, but make sure
+   if (sequence.isRing())
+   {
+      try
+      {
+         areaContainsPoint =
+            geos::algorithm::PointLocation::isInRing(zero, &sequence);
+      }
+      catch (const std::exception&)
+      {
+         logger_->trace("Invalid area sequence");
+      }
+   }
+
+   return areaContainsPoint;
+}
+
 units::angle::degrees<double>
 GetAngle(double lat1, double lon1, double lat2, double lon2)
 {

scwx-qt/source/scwx/qt/util/geographic_lib.hpp

@@ -1,5 +1,9 @@
 #pragma once
 
+#include <scwx/common/geographic.hpp>
+
+#include <vector>
+
 #include <GeographicLib/Geodesic.hpp>
 #include <units/angle.h>
 #include <units/length.h>
@@ -20,6 +24,18 @@ namespace GeographicLib
  */
 const ::GeographicLib::Geodesic& DefaultGeodesic();
 
+/**
+ * Determine if an area/ring, oriented in either direction, contains a point. A
+ * point lying on the area boundary is considered to be inside the area.
+ *
+ * @param [in] area A vector of Coordinates representing the area
+ * @param [in] point The point to check against the area
+ *
+ * @return true if point is inside the area
+ */
+bool AreaContainsPoint(const std::vector<common::Coordinate>& area,
+                       const common::Coordinate&              point);
+
 /**
  * Get the angle between two points.
  *