mirror of
https://github.com/ciphervance/supercell-wx.git
synced 2025-10-30 22:40:05 +00:00
315 lines
9.5 KiB
C++
315 lines
9.5 KiB
C++
#include <scwx/qt/util/geographic_lib.hpp>
|
|
#include <scwx/util/logger.hpp>
|
|
|
|
#include <numbers>
|
|
#include <cmath>
|
|
|
|
#include <GeographicLib/Gnomonic.hpp>
|
|
#include <geos/algorithm/PointLocation.h>
|
|
#include <geos/operation/distance/DistanceOp.h>
|
|
#include <geos/geom/CoordinateSequence.h>
|
|
#include <geos/geom/GeometryFactory.h>
|
|
|
|
namespace scwx
|
|
{
|
|
namespace qt
|
|
{
|
|
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_ {
|
|
::GeographicLib::Constants::WGS84_a(),
|
|
::GeographicLib::Constants::WGS84_f()};
|
|
|
|
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& ex)
|
|
{
|
|
logger_->warn("Invalid area sequence. {}", ex.what());
|
|
}
|
|
}
|
|
|
|
return areaContainsPoint;
|
|
}
|
|
|
|
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)
|
|
{
|
|
double azi1;
|
|
double azi2;
|
|
DefaultGeodesic().Inverse(lat1, lon1, lat2, lon2, azi1, azi2);
|
|
|
|
return units::angle::degrees<double> {azi1};
|
|
}
|
|
|
|
common::Coordinate GetCoordinate(const common::Coordinate& center,
|
|
units::angle::degrees<double> angle,
|
|
units::length::meters<double> distance)
|
|
{
|
|
double latitude;
|
|
double longitude;
|
|
|
|
DefaultGeodesic().Direct(center.latitude_,
|
|
center.longitude_,
|
|
angle.value(),
|
|
distance.value(),
|
|
latitude,
|
|
longitude);
|
|
|
|
return {latitude, longitude};
|
|
}
|
|
|
|
common::Coordinate GetCoordinate(const common::Coordinate& center,
|
|
units::meters<double> i,
|
|
units::meters<double> j)
|
|
{
|
|
// Calculate polar coordinates based on i and j
|
|
const double angle =
|
|
std::atan2(i.value(), j.value()) * 180.0 / std::numbers::pi;
|
|
const double range =
|
|
std::sqrt(i.value() * i.value() + j.value() * j.value());
|
|
|
|
double latitude;
|
|
double longitude;
|
|
|
|
DefaultGeodesic().Direct(
|
|
center.latitude_, center.longitude_, angle, range, latitude, longitude);
|
|
|
|
return {latitude, longitude};
|
|
}
|
|
|
|
units::length::meters<double>
|
|
GetDistance(double lat1, double lon1, double lat2, double lon2)
|
|
{
|
|
double distance;
|
|
DefaultGeodesic().Inverse(lat1, lon1, lat2, lon2, distance);
|
|
|
|
return units::length::meters<double> {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<double>
|
|
GetDistanceAreaPoint(const std::vector<common::Coordinate>& area,
|
|
const common::Coordinate& point)
|
|
{
|
|
// 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)
|
|
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);
|
|
// 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))
|
|
{
|
|
useCentroid = true;
|
|
}
|
|
sequence.add(x, y);
|
|
}
|
|
|
|
units::length::meters<double> distance;
|
|
|
|
if (useCentroid)
|
|
{
|
|
common::Coordinate centroid = common::GetCentroid(area);
|
|
distance = GetDistance(point.latitude_,
|
|
point.longitude_,
|
|
centroid.latitude_,
|
|
centroid.longitude_);
|
|
}
|
|
else if (GnomonicAreaContainsCenter(sequence))
|
|
{
|
|
distance = units::length::meters<double>(0);
|
|
}
|
|
else
|
|
{
|
|
// Get the closes point on the geometry
|
|
auto geometryFactory = geos::geom::GeometryFactory::getDefaultInstance();
|
|
auto lineString = geometryFactory->createLineString(sequence);
|
|
auto zeroPoint = geometryFactory->createPoint(zero);
|
|
|
|
std::unique_ptr<geos::geom::CoordinateSequence> closestPoints =
|
|
geos::operation::distance::DistanceOp::nearestPoints(lineString.get(),
|
|
zeroPoint.get());
|
|
|
|
double closestLat;
|
|
double closestLon;
|
|
|
|
gnomonic.Reverse(point.latitude_,
|
|
point.longitude_,
|
|
closestPoints->getX(0),
|
|
closestPoints->getY(0),
|
|
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;
|
|
}
|
|
|
|
units::length::meters<double>
|
|
GetRadarBeamAltititude(units::length::meters<double> range,
|
|
units::angle::degrees<double> elevation,
|
|
units::length::meters<double> height)
|
|
{
|
|
static const units::length::meters<double> earthRadius {6367444 * 4 / 3};
|
|
|
|
height += earthRadius;
|
|
|
|
const double elevationRadians =
|
|
units::angle::radians<double>(elevation).value();
|
|
const auto altitudeSquared =
|
|
(range * range + height * height +
|
|
2 * range * height * std::sin(elevationRadians));
|
|
|
|
return units::length::meters<double>(std::sqrt(altitudeSquared.value())) -
|
|
earthRadius;
|
|
}
|
|
|
|
} // namespace GeographicLib
|
|
} // namespace util
|
|
} // namespace qt
|
|
} // namespace scwx
|