Changed code to find a distance, useful for both in range and general distance.

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

@@ -31,6 +31,40 @@ const ::GeographicLib::Geodesic& DefaultGeodesic()
    return geodesic_;
 }
 
+bool GnomonicAreaContainsCenter(geos::geom::CoordinateSequence sequence)
+{
+
+   // Cannot have an area with just two points
+   if (sequence.size() <= 2 ||
+       (sequence.size() == 3 && sequence.front() == sequence.back()))
+   {
+      return false;
+   }
+   bool areaContainsPoint = false;
+   geos::geom::CoordinateXY zero {};
+   // 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;
+}
+
 bool AreaContainsPoint(const std::vector<common::Coordinate>& area,
                        const common::Coordinate&              point)
 {
@@ -146,60 +180,42 @@ GetDistance(double lat1, double lon1, double lat2, double lon2)
    return units::length::meters<double> {distance};
 }
 
-bool AreaInRangeOfPoint(const std::vector<common::Coordinate>& area,
+/*
-                        const common::Coordinate&              point,
+ * Uses the gnomonic projection to determine if the area is in the radius.
-                        const units::length::meters<double>    distance)
+ *
-{
+ * The basic algorithm is as follows:
-   /*
+ *    - Get a gnomonic projection centered on the point of the area
-   Uses the gnomonic projection to determine if the area is in the radius.
+ *    - Find the point on the area which is closest to the center
-
+ *    - Convert the closest point back to latitude and longitude.
-   The first property needed to make this work is that great circles become
+ *    - Find the distance form the closest point to the point.
-   lines in the projection.
+ *
-   The other key property needed to make this work is described bellow
+ * The first property needed to make this work is that great circles become
-      R1 and R2 are the distances from the center point to two points
+ * lines in the projection, which allows the area to be converted to strait
-      on the (non-flat) Earth.
+ * lines. This is generally true for gnomic projections.
-      R1' and R2' are the distances from the center point to the same
+ *
-      two points in the gnomonic projection.
+ * The second property needed to make this work is that a point further away
-      if R1 > R2 then
+ * on the geodesic must be further away on the projection. This means that
-         R1' > R2'
+ * the closes point on the projection is also the closest point on the geodesic.
-      else if R1 < R2 then
+ * This holds for spherical geodesics and is an approximation non spherical
-         R1' < R2'
+ * geodesics.
-      else if R1 == R2 then
+ *
-         R1' == R2'
+ * This algorithm only works if the area is fully on the hemisphere centered
-
+ * on the point. Otherwise, this falls back to centroid based distances.
-      This can also be written as:
+ *
-      r(d) is a function that takes the distance on Earth and converts it to a
+ * If the point is inside the area, 0 is always returned.
-      distance on the projection.
-      R1' = r(R1), R2' = r(R2)
-      r(d) is increasing
-
-   In this case, R1 is a point the radius away from the center, and R2 is a
-   (all of the) point(s) on the edge of the area. This means that if the edge
-   is in the radius R1' on the projection, it is in the radius R1 on the Earth.
-
-   On a spherical geodesic this works fine. R is the radius of Earth. We are
-   also only concerned with points less than a hemisphere away, therefore
-   0 < R1,R2 < pi/2 * R (quarter of circumference because the point is in the
-   center of the hemisphere)
-      r(d) = R * tan(d / R) {0 < d < pi/2 * R}
-      tan(d / R) is increasing for {0 < d < pi/2 * R}
-
-   On non spherical geodesics, this may not work perfectly, but should be a
-   close approximation.
  */
-   // Cannot have an area with just two points
+units::length::meters<double>
-   if (area.size() <= 2 || (area.size() == 3 && area.front() == area.back()))
+GetDistanceAreaPoint(const std::vector<common::Coordinate>& area,
-   {
+                     const common::Coordinate&              point)
-      return false;
+{
-   }
+   // Ensure that the same geodesic is used here as is for the distance
-
-   // Ensure that the same geodesic is used here as is for the radius
    // calculation
    ::GeographicLib::Gnomonic gnomonic =
       ::GeographicLib::Gnomonic(DefaultGeodesic());
    geos::geom::CoordinateSequence sequence {};
    double                         x;
    double                         y;
+   bool                           useCentroid = false;
 
    // Using a gnomonic projection with the test point as the center
    // latitude/longitude, the projected test point will be at (0, 0)
@@ -214,73 +230,64 @@ bool AreaInRangeOfPoint(const std::vector<common::Coordinate>& area,
                        areaCoordinate.longitude_,
                        x,
                        y);
-      // Check if the current point is the hemisphere centered on the point
+      // Check if the current point is in the hemisphere centered on the point
+      // if not, fall back to using centroid.
       if (std::isnan(x) || std::isnan(y))
       {
-         return false;
+         useCentroid = true;
       }
       sequence.add(x, y);
    }
 
-   // get a point on the circle with the radius of the range in lat lon.
+   units::length::meters<double> distance;
-   // Has the point be in the general direction of the area, which may help with
+
-   // non spherical geodesics
+   if (useCentroid)
-   units::angle::degrees<double> angle = GetAngle(
+   {
-      point.latitude_, point.longitude_, area[0].latitude_, area[0].longitude_);
+      common::Coordinate centroid = common::GetCentroid(area);
-   common::Coordinate radiusPoint = GetCoordinate(point, angle, distance);
+      distance = GetDistance(point.latitude_,
-   // get the radius in gnomonic projection
-   gnomonic.Forward(point.latitude_,
                              point.longitude_,
-                    radiusPoint.latitude_,
+                             centroid.latitude_,
-                    radiusPoint.longitude_,
+                             centroid.longitude_);
-                    x,
-                    y);
-   // radius is greater than quarter circumference of the Earth, but the area
-   // is closer, so it is in range.
-   if (std::isnan(x) || std::isnan(y))
-   {
-      return true;
    }
-   double gnomonicRadius = std::sqrt(x * x + y * y);
+   else if (GnomonicAreaContainsCenter(sequence))
-
-   // If the sequence is not a ring, add the first point again for closure
-   if (!sequence.isRing())
    {
-      sequence.add(sequence.front(), false);
+      distance = units::length::meters<double>(0);
    }
-
+   else
-   // The sequence should be a ring at this point, but make sure
-   if (sequence.isRing())
    {
-      try
+      // Get the closes point on the geometry
-      {
+      auto geometryFactory = geos::geom::GeometryFactory::getDefaultInstance();
-         if (geos::algorithm::PointLocation::isInRing(zero, &sequence))
+      auto lineString      = geometryFactory->createLineString(sequence);
-         {
-            return true;
-         }
-         else if (distance > units::length::meters<double>(0))
-         {
-            // Calculate the distance the area is from the point via conversion
-            // to a polygon.
-            auto geometryFactory =
-               geos::geom::GeometryFactory::getDefaultInstance();
-            auto linearRing = geometryFactory->createLinearRing(sequence);
-            auto polygon =
-               geometryFactory->createPolygon(std::move(linearRing));
 
       geos::algorithm::distance::PointPairDistance distancePair;
       geos::algorithm::distance::DistanceToPoint::computeDistance(
-               *polygon, zero, distancePair);
+         *lineString, zero, distancePair);
-            return gnomonicRadius >= distancePair.getDistance();
-         }
-      }
-      catch (const std::exception&)
-      {
-         logger_->trace("Invalid area sequence");
-      }
-   }
 
-   return false;
+      geos::geom::CoordinateXY closestPoint = distancePair.getCoordinate(0);
+
+      double closestLat;
+      double closestLon;
+
+      gnomonic.Reverse(point.latitude_,
+                       point.longitude_,
+                       closestPoint.x,
+                       closestPoint.y,
+                       closestLat,
+                       closestLon);
+
+      distance = GetDistance(point.latitude_,
+                             point.longitude_,
+                             closestLat,
+                             closestLon);
+   }
+   return distance;
+}
+
+bool AreaInRangeOfPoint(const std::vector<common::Coordinate>& area,
+                        const common::Coordinate&              point,
+                        const units::length::meters<double>    distance)
+{
+    return GetDistanceAreaPoint(area, point) <= distance;
 }
 
 } // namespace GeographicLib

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

@@ -90,14 +90,26 @@ common::Coordinate GetCoordinate(const common::Coordinate& center,
 units::length::meters<double>
 GetDistance(double lat1, double lon1, double lat2, double lon2);
 
+/**
+ * Get the distance from an area to a point. If the area is less than a quarter
+ * radius of the Earth away, this is the closest distance between the area and
+ * the point. Otherwise it is the distance from the centroid of the area to the
+ * point. Finally, if the point is in the area, it is always 0.
+ *
+ * @param [in] area A vector of Coordinates representing the area
+ * @param [in] point The point to check against the area
+ *
+ * @return true if area is inside the radius of the point
+ */
+units::length::meters<double>
+GetDistanceAreaPoint(const std::vector<common::Coordinate>& area,
+                     const common::Coordinate&              point);
 
 /**
  * Determine if an area/ring, oriented in either direction, is within a
  * distance of a point. A point lying on the area boundary is considered to be
  * inside the area, and thus always in range. Any part of the area being inside
- * the radius counts as inside.
+ * the radius counts as inside. Uses GetDistanceAreaPoint to get the distance.
- * This is limited to having the area be in the same hemisphere centered on
- * the point, and radices up to a quarter of the circumference of the Earth.
  *
  * @param [in] area A vector of Coordinates representing the area
  * @param [in] point The point to check against the area