diff --git a/scwx-qt/source/scwx/qt/util/geographic_lib.cpp b/scwx-qt/source/scwx/qt/util/geographic_lib.cpp index 602f03af..2557fd15 100644 --- a/scwx-qt/source/scwx/qt/util/geographic_lib.cpp +++ b/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& area, const common::Coordinate& point) { @@ -146,60 +180,42 @@ GetDistance(double lat1, double lon1, double lat2, double lon2) return units::length::meters {distance}; } -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 basic algorithm is as follows: + * - Get a gnomonic projection centered on the point of the area + * - Find the point on the area which is closest to the center + * - Convert the closest point back to latitude and longitude. + * - Find the distance form the closest point to the point. + * + * The first property needed to make this work is that great circles become + * lines in the projection, which allows the area to be converted to strait + * lines. This is generally true for gnomic projections. + * + * The second property needed to make this work is that a point further away + * on the geodesic must be further away on the projection. This means that + * the closes point on the projection is also the closest point on the geodesic. + * This holds for spherical geodesics and is an approximation non spherical + * geodesics. + * + * This algorithm only works if the area is fully on the hemisphere centered + * on the point. Otherwise, this falls back to centroid based distances. + * + * If the point is inside the area, 0 is always returned. + */ +units::length::meters +GetDistanceAreaPoint(const std::vector& area, + const common::Coordinate& point) { - /* - 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; - } - - // Ensure that the same geodesic is used here as is for the radius + // Ensure that the same geodesic is used here as is for the distance // 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& 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. - // 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_, - point.longitude_, - radiusPoint.latitude_, - radiusPoint.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)) + units::length::meters distance; + + if (useCentroid) { - return true; + common::Coordinate centroid = common::GetCentroid(area); + distance = GetDistance(point.latitude_, + point.longitude_, + centroid.latitude_, + centroid.longitude_); } - 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()) + else if (GnomonicAreaContainsCenter(sequence)) { - sequence.add(sequence.front(), false); + distance = units::length::meters(0); } - - // The sequence should be a ring at this point, but make sure - if (sequence.isRing()) + else { - try - { - if (geos::algorithm::PointLocation::isInRing(zero, &sequence)) - { - return true; - } - else if (distance > units::length::meters(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)); + // Get the closes point on the geometry + auto geometryFactory = geos::geom::GeometryFactory::getDefaultInstance(); + auto lineString = geometryFactory->createLineString(sequence); - geos::algorithm::distance::PointPairDistance distancePair; - geos::algorithm::distance::DistanceToPoint::computeDistance( - *polygon, zero, distancePair); - return gnomonicRadius >= distancePair.getDistance(); - } - } - catch (const std::exception&) - { - logger_->trace("Invalid area sequence"); - } + geos::algorithm::distance::PointPairDistance distancePair; + geos::algorithm::distance::DistanceToPoint::computeDistance( + *lineString, zero, distancePair); + + 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; +} - return false; +bool AreaInRangeOfPoint(const std::vector& area, + const common::Coordinate& point, + const units::length::meters distance) +{ + return GetDistanceAreaPoint(area, point) <= distance; } } // namespace GeographicLib diff --git a/scwx-qt/source/scwx/qt/util/geographic_lib.hpp b/scwx-qt/source/scwx/qt/util/geographic_lib.hpp index a7509442..5038d9a9 100644 --- a/scwx-qt/source/scwx/qt/util/geographic_lib.hpp +++ b/scwx-qt/source/scwx/qt/util/geographic_lib.hpp @@ -90,14 +90,26 @@ common::Coordinate GetCoordinate(const common::Coordinate& center, units::length::meters 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 +GetDistanceAreaPoint(const std::vector& 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. - * 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. + * the radius counts as inside. Uses GetDistanceAreaPoint to get the distance. * * @param [in] area A vector of Coordinates representing the area * @param [in] point The point to check against the area