From 2bc971eb94a4ba70d389d245dc35f326aeb0326b Mon Sep 17 00:00:00 2001 From: AdenKoperczak Date: Mon, 1 Jul 2024 16:16:54 -0400 Subject: [PATCH] added comments and improved code layout --- .../source/scwx/qt/util/geographic_lib.cpp | 91 +++++++++++++++---- .../source/scwx/qt/util/geographic_lib.hpp | 2 + 2 files changed, 73 insertions(+), 20 deletions(-) diff --git a/scwx-qt/source/scwx/qt/util/geographic_lib.cpp b/scwx-qt/source/scwx/qt/util/geographic_lib.cpp index bc4066ce..602f03af 100644 --- a/scwx-qt/source/scwx/qt/util/geographic_lib.cpp +++ b/scwx-qt/source/scwx/qt/util/geographic_lib.cpp @@ -14,7 +14,8 @@ namespace scwx { namespace qt { -namespace util { +namespace util +{ namespace GeographicLib { @@ -149,15 +150,52 @@ bool AreaInRangeOfPoint(const std::vector& area, const common::Coordinate& point, const units::length::meters distance) { + /* + Uses the gnomonic projection to determine if the area is in the radius. + + The first property needed to make this work is that great circles become + lines in the projection. + The other key property needed to make this work is described bellow + R1 and R2 are the distances from the center point to two points + on the (non-flat) Earth. + R1' and R2' are the distances from the center point to the same + two points in the gnomonic projection. + if R1 > R2 then + R1' > R2' + else if R1 < R2 then + R1' < R2' + else if R1 == R2 then + R1' == R2' + + This can also be written as: + r(d) is a function that takes the distance on Earth and converts it to a + 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 if (area.size() <= 2 || (area.size() == 3 && area.front() == area.back())) { return false; } - - - ::GeographicLib::Gnomonic gnomonic = + // 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; @@ -176,11 +214,19 @@ bool AreaInRangeOfPoint(const std::vector& area, areaCoordinate.longitude_, x, y); + // Check if the current point is the hemisphere centered on the point + if (std::isnan(x) || std::isnan(y)) + { + return false; + } sequence.add(x, y); } // get a point on the circle with the radius of the range in lat lon. - units::angle::degrees angle = units::angle::degrees(0); + // Has the point be in the general direction of the area, which may help with + // non spherical geodesics + units::angle::degrees angle = GetAngle( + point.latitude_, point.longitude_, area[0].latitude_, area[0].longitude_); common::Coordinate radiusPoint = GetCoordinate(point, angle, distance); // get the radius in gnomonic projection gnomonic.Forward(point.latitude_, @@ -189,7 +235,13 @@ bool AreaInRangeOfPoint(const std::vector& area, radiusPoint.longitude_, x, y); - double gnomonicRadius = sqrt(x * x + y * 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); // If the sequence is not a ring, add the first point again for closure if (!sequence.isRing()) @@ -206,22 +258,21 @@ bool AreaInRangeOfPoint(const std::vector& area, { return true; } - - // Calculate the distance the point is from the output - geos::algorithm::distance::PointPairDistance distancePair; - auto geometryFactory = - geos::geom::GeometryFactory::getDefaultInstance(); - auto linearRing = geometryFactory->createLinearRing(sequence); - auto polygon = - geometryFactory->createPolygon(std::move(linearRing)); - geos::algorithm::distance::DistanceToPoint::computeDistance(*polygon, - zero, - distancePair); - if (gnomonicRadius > distancePair.getDistance()) + else if (distance > units::length::meters(0)) { - return true; - } + // 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); + return gnomonicRadius >= distancePair.getDistance(); + } } catch (const std::exception&) { diff --git a/scwx-qt/source/scwx/qt/util/geographic_lib.hpp b/scwx-qt/source/scwx/qt/util/geographic_lib.hpp index 93469802..a7509442 100644 --- a/scwx-qt/source/scwx/qt/util/geographic_lib.hpp +++ b/scwx-qt/source/scwx/qt/util/geographic_lib.hpp @@ -96,6 +96,8 @@ GetDistance(double lat1, double lon1, double lat2, double lon2); * 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. + * 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