XAML-Map-Control/MapControl/Shared/MapProjection.cs

using System;
#if WPF
using System.Windows;
using System.Windows.Media;
#elif AVALONIA
using Avalonia;
#endif

namespace MapControl
{
    /// <summary>
    /// Implements a map projection, a transformation between geographic coordinates,
    /// i.e. latitude and longitude in degrees, and cartesian map coordinates in meters.
    /// </summary>
#if UWP || WINUI
    [Windows.Foundation.Metadata.CreateFromString(MethodName = "Parse")]
#else
    [System.ComponentModel.TypeConverter(typeof(MapProjectionConverter))]
#endif
    public abstract class MapProjection
    {
        public const double Wgs84EquatorialRadius = 6378137d;
        public const double Wgs84Flattening = 1d / 298.257223563;
        public const double Wgs84MeterPerDegree = Wgs84EquatorialRadius * Math.PI / 180d;

        public static MapProjectionFactory Factory
        {
            get => field ??= new MapProjectionFactory();
            set;
        }

        /// <summary>
        /// Creates a MapProjection instance from a CRS identifier string.
        /// </summary>
        public static MapProjection Parse(string crsId)
        {
            return Factory.GetProjection(crsId);
        }

        /// <summary>
        /// Gets the WMS 1.3.0 CRS identifier.
        /// </summary>
        public string CrsId { get; protected set; } = "";

        /// <summary>
        /// Indicates whether the projection is normal cylindrical, see
        /// https://en.wikipedia.org/wiki/Map_projection#Normal_cylindrical.
        /// </summary>
        public bool IsNormalCylindrical { get; protected set; }

        /// <summary>
        /// The earth ellipsoid semi-major axis, or spherical earth radius respectively, in meters.
        /// </summary>
        public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius;

        public double MeterPerDegree => EquatorialRadius * Math.PI / 180d;

        private Location center;
        private bool updateCenter;

        /// <summary>
        /// Gets or sets an optional projection center. If the property is set to a non-null value,
        /// it overrides the projection center set by MapBase.Center or MapBase.ProjectionCenter.
        /// </summary>
        public Location Center
        {
            get => center ??= new Location();
            set
            {
                updateCenter = true;

                if (value != null)
                {
                    var longitude = Location.NormalizeLongitude(value.Longitude);

                    SetCenter(value.LongitudeEquals(longitude) ? value : new Location(value.Latitude, longitude));

                    updateCenter = false;
                }
            }
        }

        protected void EnableCenterUpdates()
        {
            updateCenter = true;
        }

        /// <summary>
        /// Called by MapBase.UpdateTransform().
        /// </summary>
        internal void SetCenter(Location value)
        {
            if (updateCenter)
            {
                if (center == null || !center.Equals(value))
                {
                    center = value;
                    CenterChanged();
                }
            }
        }

        protected virtual void CenterChanged() { }

        /// <summary>
        /// Gets the relative scale at the specified geographic coordinates.
        /// The returned Matrix represents the local distortion of the map projection.
        /// </summary>
        public abstract Matrix RelativeScale(double latitude, double longitude);

        /// <summary>
        /// Transforms geographic coordinates to a Point in projected map coordinates.
        /// Returns null when the location can not be transformed.
        /// </summary>
        public abstract Point? LocationToMap(double latitude, double longitude);

        /// <summary>
        /// Transforms projected map coordinates to a Location in geographic coordinates.
        /// Returns null when the coordinates can not be transformed.
        /// </summary>
        public abstract Location MapToLocation(double x, double y);

        /// <summary>
        /// Gets the relative map scale at the specified geographic Location.
        /// </summary>
        public Matrix RelativeScale(Location location) => RelativeScale(location.Latitude, location.Longitude);

        /// <summary>
        /// Transforms a Location in geographic coordinates to a Point in projected map coordinates.
        /// Returns null when the Location can not be transformed.
        /// </summary>
        public Point? LocationToMap(Location location) => LocationToMap(location.Latitude, location.Longitude);

        /// <summary>
        /// Transforms a Point in projected map coordinates to a Location in geographic coordinates.
        /// Returns null when the Point can not be transformed.
        /// </summary>
        public Location MapToLocation(Point point) => MapToLocation(point.X, point.Y);

        /// <summary>
        /// Transforms a BoundingBox in geographic coordinates to a Rect in projected map coordinates
        /// with an optional transform Matrix when the BoundingBox is not projected axis-aligned.
        /// Returns (null, null) when the BoundingBox can not be transformed.
        /// </summary>
        public (Rect?, Matrix?) BoundingBoxToMap(BoundingBox boundingBox)
        {
            Rect? rect = null;
            Matrix? transform = null;
            var sw = LocationToMap(boundingBox.South, boundingBox.West);
            var ne = LocationToMap(boundingBox.North, boundingBox.East);

            if (sw.HasValue && ne.HasValue)
            {
                if (boundingBox.ProjectAxisAligned)
                {
                    rect = new Rect(sw.Value, ne.Value);
                }
                else
                {
                    var se = LocationToMap(boundingBox.South, boundingBox.East);
                    var nw = LocationToMap(boundingBox.North, boundingBox.West);

                    if (se.HasValue && nw.HasValue)
                    {
                        var south = new Point((sw.Value.X + se.Value.X) / 2d, (sw.Value.Y + se.Value.Y) / 2d); // south midpoint
                        var north = new Point((nw.Value.X + ne.Value.X) / 2d, (nw.Value.Y + ne.Value.Y) / 2d); // north midpoint
                        var west = new Point((nw.Value.X + sw.Value.X) / 2d, (nw.Value.Y + sw.Value.Y) / 2d);  // west midpoint
                        var east = new Point((ne.Value.X + se.Value.X) / 2d, (ne.Value.Y + se.Value.Y) / 2d);  // east midpoint
                        var center = new Point((west.X + east.X) / 2d, (west.Y + east.Y) / 2d); // midpoint of segment west-east
                        var dx1 = east.X - west.X;
                        var dy1 = east.Y - west.Y;
                        var dx2 = north.X - south.X;
                        var dy2 = north.Y - south.Y;
                        var width = Math.Sqrt(dx1 * dx1 + dy1 * dy1); // distance west-east
                        var height = Math.Sqrt(dx2 * dx2 + dy2 * dy2); // distance south-north

                        rect = new Rect(center.X - width / 2d, center.Y - height / 2d, width, height);

                        if (dy1 != 0d || dx2 != 0d)
                        {
                            // Skew matrix with skewX = Atan(-dx2 / dy2) and skewY = Atan(-dy1 / dx1).
                            //
                            transform = new Matrix(1d, -dy1 / dx1, -dx2 / dy2, 1d, 0d, 0d);
                        }
                    }
                }
            }

            return (rect, transform);
        }

        /// <summary>
        /// Transforms a Rect in projected map coordinates to a BoundingBox in geographic coordinates.
        /// Returns null when the Rect can not be transformed.
        /// </summary>
        public BoundingBox MapToBoundingBox(Rect rect, bool axisAligned = false)
        {
            var sw = MapToLocation(rect.X, rect.Y);
            var ne = MapToLocation(rect.X + rect.Width, rect.Y + rect.Height);

            return sw != null && ne != null
                ? new BoundingBox(sw.Latitude, sw.Longitude, ne.Latitude, ne.Longitude, axisAligned)
                : null;
        }

        public override string ToString()
        {
            return CrsId;
        }
    }
}