annotations/Azaion.Common/Extensions/GeoUtils.cs

using Azaion.Common.DTO;

namespace Azaion.Common.Extensions;

public static class GeoUtils
{
    private const double EARTH_RADIUS = 6378137;

    public static (int x, int y) WorldToTilePos(double lat, double lon, int zoom)
    {
        var latRad = lat * Math.PI / 180.0;
        var n = Math.Pow(2.0, zoom);
        var xTile = (int)Math.Floor((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 = EARTH_RADIUS * 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 / EARTH_RADIUS;
        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 / EARTH_RADIUS) * (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 / (EARTH_RADIUS * 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);
    }
}