satellite-provider/SatelliteProvider.Common/Utils/GeoUtils.cs

using SatelliteProvider.Common.DTO;

namespace SatelliteProvider.Common.Utils;

public static class GeoUtils
{
    public const double EarthRadiusMeters = 6378137d;
    public const double EarthEquatorialCircumferenceMeters = 40075016.686d;
    public const double MetersPerDegreeLatitude = 111000d;

    public static (int x, int y) WorldToTilePos(GeoPoint point, int zoom)
    {
        var latRad = point.Lat * Math.PI / 180.0;
        var n = Math.Pow(2.0, zoom);
        var xTile = (int)Math.Floor((point.Lon + 180.0) / 360.0 * n);
        var yTile = (int)Math.Floor((1.0 - Math.Log(Math.Tan(latRad) + 1.0 / Math.Cos(latRad)) / Math.PI) / 2.0 * n);
        return (xTile, yTile);
    }

    public static double ToRadians(double degrees) => degrees * Math.PI / 180.0;
    public static double ToDegrees(double radians) => radians * 180.0 / Math.PI;

    public static Direction DirectionTo(this GeoPoint p1, GeoPoint p2)
    {
        var lat1Rad = ToRadians(p1.Lat);
        var lat2Rad = ToRadians(p2.Lat);
        var dLon = ToRadians(p2.Lon - p1.Lon);
        var dLat = ToRadians(p2.Lat - p1.Lat);

        var a = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) +
                Math.Cos(lat1Rad) * Math.Cos(lat2Rad) *
                Math.Sin(dLon / 2) * Math.Sin(dLon / 2);
        var c = 2 * Math.Asin(Math.Sqrt(a));
        var distance = EarthRadiusMeters * c;

        var y = Math.Sin(dLon) * Math.Cos(lat2Rad);
        var x = Math.Cos(lat1Rad) * Math.Sin(lat2Rad) -
                   Math.Sin(lat1Rad) * Math.Cos(lat2Rad) * Math.Cos(dLon);
        var azimuthRadians = Math.Atan2(y, x);
        var azimuth = (ToDegrees(azimuthRadians) + 360) % 360;

        return new Direction
        {
            Distance = distance,
            Azimuth = azimuth
        };
    }

    public static GeoPoint GoDirection(this GeoPoint startPoint, Direction direction)
    {
        var angularDistance = direction.Distance / EarthRadiusMeters;
        var azimuthRadians = ToRadians(direction.Azimuth);
        var startLatRad = ToRadians(startPoint.Lat);
        var startLonRad = ToRadians(startPoint.Lon);

        var destLatRad = Math.Asin(Math.Sin(startLatRad) * Math.Cos(angularDistance) +
                                   Math.Cos(startLatRad) * Math.Sin(angularDistance) * Math.Cos(azimuthRadians));

        var destLonRad = startLonRad + Math.Atan2(Math.Sin(azimuthRadians) * Math.Sin(angularDistance) * Math.Cos(startLatRad),
            Math.Cos(angularDistance) - Math.Sin(startLatRad) * Math.Sin(destLatRad));

        return new GeoPoint(ToDegrees(destLatRad), ToDegrees(destLonRad));
    }

    public static GeoPoint TileToWorldPos(int x, int y, int zoom)
    {
        var n = Math.Pow(2.0, zoom);
        var lonDeg = x / n * 360.0 - 180.0;
        var latRad = Math.Atan(Math.Sinh(Math.PI * (1.0 - 2.0 * y / n)));
        var latDeg = latRad * 180.0 / Math.PI;
        return new GeoPoint(latDeg, lonDeg);
    }

    public static (double minLat, double maxLat, double minLon, double maxLon) GetBoundingBox(GeoPoint centerGeoPoint, double radiusM)
    {
        var latRad = centerGeoPoint.Lat * Math.PI / 180.0;

        var latDiff = (radiusM / EarthRadiusMeters) * (180.0 / Math.PI);
        var minLat = Math.Max(centerGeoPoint.Lat - latDiff, -90.0);
        var maxLat = Math.Min(centerGeoPoint.Lat + latDiff, 90.0);

        var lonDiff = (radiusM / (EarthRadiusMeters * Math.Cos(latRad))) * (180.0 / Math.PI);
        var minLon = Math.Max(centerGeoPoint.Lon - lonDiff, -180.0);
        var maxLon = Math.Min(centerGeoPoint.Lon + lonDiff, 180.0);

        return (minLat, maxLat, minLon, maxLon);
    }

    public static List<GeoPoint> CalculateIntermediatePoints(GeoPoint start, GeoPoint end, double maxSpacingMeters)
    {
        var direction = start.DirectionTo(end);
        var distance = direction.Distance;

        if (distance <= maxSpacingMeters)
        {
            return new List<GeoPoint>();
        }

        var numSegments = (int)Math.Ceiling(distance / maxSpacingMeters);
        var actualSpacing = distance / numSegments;

        var intermediatePoints = new List<GeoPoint>();

        for (int i = 1; i < numSegments; i++)
        {
            var segmentDistance = actualSpacing * i;
            var intermediateDirection = new Direction
            {
                Distance = segmentDistance,
                Azimuth = direction.Azimuth
            };
            var intermediatePoint = start.GoDirection(intermediateDirection);
            intermediatePoints.Add(intermediatePoint);
        }

        return intermediatePoints;
    }

    public static double CalculateDistance(GeoPoint p1, GeoPoint p2)
    {
        return p1.DirectionTo(p2).Distance;
    }

    public static GeoPoint CalculateCenter(GeoPoint northWest, GeoPoint southEast)
    {
        var centerLat = (northWest.Lat + southEast.Lat) / 2.0;
        var centerLon = (northWest.Lon + southEast.Lon) / 2.0;
        return new GeoPoint(centerLat, centerLon);
    }
}