diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs index 52016c77..3d5911bc 100644 --- a/MapControl/Shared/AzimuthalEquidistantProjection.cs +++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs @@ -10,7 +10,7 @@ using System.Windows; namespace MapControl { /// - /// Transforms map coordinates according to the Azimuthal Equidistant Projection. + /// Spherical Azimuthal Equidistant Projection. /// public class AzimuthalEquidistantProjection : AzimuthalProjection { @@ -35,7 +35,7 @@ namespace MapControl GetAzimuthDistance(ProjectionCenter, location, out azimuth, out distance); - distance *= Wgs84EquatorialRadius; + distance *= TrueScale * 180d / Math.PI; return new Point(distance * Math.Sin(azimuth), distance * Math.Cos(azimuth)); } @@ -48,7 +48,7 @@ namespace MapControl } var azimuth = Math.Atan2(point.X, point.Y); - var distance = Math.Sqrt(point.X * point.X + point.Y * point.Y) / Wgs84EquatorialRadius; + var distance = Math.Sqrt(point.X * point.X + point.Y * point.Y) / (TrueScale * 180d / Math.PI); return GetLocation(ProjectionCenter, azimuth, distance); } diff --git a/MapControl/Shared/AzimuthalProjection.cs b/MapControl/Shared/AzimuthalProjection.cs index 0153ee77..8a8a865f 100644 --- a/MapControl/Shared/AzimuthalProjection.cs +++ b/MapControl/Shared/AzimuthalProjection.cs @@ -67,7 +67,7 @@ namespace MapControl } /// - /// Calculates azimuth and distance in radians from location1 to location2. + /// Calculates azimuth and spherical distance in radians from location1 to location2. /// The returned distance has to be multiplied with an appropriate earth radius. /// public static void GetAzimuthDistance(Location location1, Location location2, out double azimuth, out double distance) @@ -89,7 +89,7 @@ namespace MapControl } /// - /// Calculates the Location of the point given by azimuth and distance in radians from location. + /// Calculates the Location of the point given by azimuth and spherical distance in radians from location. /// public static Location GetLocation(Location location, double azimuth, double distance) { diff --git a/MapControl/Shared/BoundingBox.cs b/MapControl/Shared/BoundingBox.cs index 61b244fa..6dcd8709 100644 --- a/MapControl/Shared/BoundingBox.cs +++ b/MapControl/Shared/BoundingBox.cs @@ -56,11 +56,6 @@ namespace MapControl get { return North - South; } } - public bool HasValidBounds - { - get { return South < North && West < East; } - } - public virtual BoundingBox Clone() { return new BoundingBox(South, West, North, East); diff --git a/MapControl/Shared/CenteredBoundingBox.cs b/MapControl/Shared/CenteredBoundingBox.cs index 36bb24a7..87722e44 100644 --- a/MapControl/Shared/CenteredBoundingBox.cs +++ b/MapControl/Shared/CenteredBoundingBox.cs @@ -2,6 +2,8 @@ // © 2018 Clemens Fischer // Licensed under the Microsoft Public License (Ms-PL) +using System; + namespace MapControl { public class CenteredBoundingBox : BoundingBox @@ -12,8 +14,8 @@ namespace MapControl public CenteredBoundingBox(Location center, double width, double height) { Center = center; - this.width = width; - this.height = height; + this.width = Math.Max(width, 0d); + this.height = Math.Max(height, 0d); } public Location Center { get; private set; } diff --git a/MapControl/Shared/EquirectangularProjection.cs b/MapControl/Shared/EquirectangularProjection.cs index 4fe1bd5e..99a746d7 100644 --- a/MapControl/Shared/EquirectangularProjection.cs +++ b/MapControl/Shared/EquirectangularProjection.cs @@ -10,7 +10,7 @@ using System.Windows; namespace MapControl { /// - /// Transforms map coordinates according to the Equirectangular Projection. + /// Equirectangular Projection. /// Longitude and Latitude values are transformed identically to X and Y. /// public class EquirectangularProjection : MapProjection @@ -23,8 +23,8 @@ namespace MapControl public EquirectangularProjection(string crsId) { CrsId = crsId; - IsCylindrical = true; - TrueScale = 1; + IsNormalCylindrical = true; + TrueScale = 1d; } public override Vector GetMapScale(Location location) diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs index dfdf9b75..31105ff7 100644 --- a/MapControl/Shared/GnomonicProjection.cs +++ b/MapControl/Shared/GnomonicProjection.cs @@ -10,7 +10,7 @@ using System.Windows; namespace MapControl { /// - /// Transforms map coordinates according to the Gnomonic Projection. + /// Spherical Gnomonic Projection. /// public class GnomonicProjection : AzimuthalProjection { @@ -36,7 +36,7 @@ namespace MapControl GetAzimuthDistance(ProjectionCenter, location, out azimuth, out distance); var mapDistance = distance < Math.PI / 2d - ? Wgs84EquatorialRadius * Math.Tan(distance) + ? Math.Tan(distance) * TrueScale * 180d / Math.PI : double.PositiveInfinity; return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth)); @@ -51,7 +51,7 @@ namespace MapControl var azimuth = Math.Atan2(point.X, point.Y); var mapDistance = Math.Sqrt(point.X * point.X + point.Y * point.Y); - var distance = Math.Atan(mapDistance / Wgs84EquatorialRadius); + var distance = Math.Atan(mapDistance / (TrueScale * 180d / Math.PI)); return GetLocation(ProjectionCenter, azimuth, distance); } diff --git a/MapControl/Shared/LocationEx.cs b/MapControl/Shared/LocationEx.cs index 1a901e28..6f45c384 100644 --- a/MapControl/Shared/LocationEx.cs +++ b/MapControl/Shared/LocationEx.cs @@ -62,7 +62,7 @@ namespace MapControl var lat2 = Math.Asin(sinLat1 * cosS12 + cosLat1 * sinS12 * cosAz1); var lon2 = lon1 + Math.Atan2(sinS12 * sinAz1, (cosLat1 * cosS12 - sinLat1 * sinS12 * cosAz1)); - return new Location(lat2 / Math.PI * 180d, lon2 / Math.PI * 180d); + return new Location(lat2 * 180d / Math.PI, lon2 * 180d / Math.PI); } public static LocationCollection CalculateMeridianLocations(this Location location, double latitude2, double resolution = 1d) diff --git a/MapControl/Shared/MapBase.cs b/MapControl/Shared/MapBase.cs index 91ec8db6..5fff98e2 100644 --- a/MapControl/Shared/MapBase.cs +++ b/MapControl/Shared/MapBase.cs @@ -348,17 +348,14 @@ namespace MapControl /// public void ZoomToBounds(BoundingBox boundingBox) { - if (boundingBox != null && boundingBox.HasValidBounds) - { - var rect = MapProjection.BoundingBoxToRect(boundingBox); - var center = new Point(rect.X + rect.Width / 2d, rect.Y + rect.Height / 2d); - var scale = Math.Min(RenderSize.Width / rect.Width, RenderSize.Height / rect.Height) - * MapProjection.TrueScale / MapProjection.PixelPerDegree; + var rect = MapProjection.BoundingBoxToRect(boundingBox); + var center = new Point(rect.X + rect.Width / 2d, rect.Y + rect.Height / 2d); + var scale = Math.Min(RenderSize.Width / rect.Width, RenderSize.Height / rect.Height) + * MapProjection.TrueScale / MapProjection.PixelPerDegree; - TargetZoomLevel = Math.Log(scale, 2d); - TargetCenter = MapProjection.PointToLocation(center); - TargetHeading = 0d; - } + TargetZoomLevel = Math.Log(scale, 2d); + TargetCenter = MapProjection.PointToLocation(center); + TargetHeading = 0d; } private void MapLayerPropertyChanged(UIElement oldLayer, UIElement newLayer) @@ -451,7 +448,7 @@ namespace MapControl if (!targetCenter.Equals(Center)) { - var targetCenterLongitude = MapProjection.IsCylindrical + var targetCenterLongitude = MapProjection.IsNormalCylindrical ? Location.NearestLongitude(targetCenter.Longitude, Center.Longitude) : targetCenter.Longitude; diff --git a/MapControl/Shared/MapImageLayer.cs b/MapControl/Shared/MapImageLayer.cs index f78637e4..e3da420c 100644 --- a/MapControl/Shared/MapImageLayer.cs +++ b/MapControl/Shared/MapImageLayer.cs @@ -256,7 +256,7 @@ namespace MapControl boundingBox = ParentMap.MapProjection.ViewportRectToBoundingBox(rect); - if (boundingBox != null && boundingBox.HasValidBounds) + if (boundingBox != null) { if (!double.IsNaN(MinLatitude) && boundingBox.South < MinLatitude) { @@ -289,7 +289,7 @@ namespace MapControl private void AdjustBoundingBox(double longitudeOffset) { - if (Math.Abs(longitudeOffset) > 180d && boundingBox != null && boundingBox.HasValidBounds) + if (Math.Abs(longitudeOffset) > 180d && boundingBox != null) { var offset = 360d * Math.Sign(longitudeOffset); @@ -300,7 +300,7 @@ namespace MapControl { var bbox = GetBoundingBox(element); - if (bbox != null && bbox.HasValidBounds) + if (bbox != null) { SetBoundingBox(element, new BoundingBox(bbox.South, bbox.West + offset, bbox.North, bbox.East + offset)); } diff --git a/MapControl/Shared/MapPanel.cs b/MapControl/Shared/MapPanel.cs index ba6cab41..762a2338 100644 --- a/MapControl/Shared/MapPanel.cs +++ b/MapControl/Shared/MapPanel.cs @@ -247,7 +247,7 @@ namespace MapControl var projection = parentMap.MapProjection; pos = projection.LocationToViewportPoint(location); - if (projection.IsCylindrical && + if (projection.IsNormalCylindrical && (pos.X < 0d || pos.X > parentMap.RenderSize.Width || pos.Y < 0d || pos.Y > parentMap.RenderSize.Height)) { @@ -311,7 +311,7 @@ namespace MapControl var center = new Point(rect.X + rect.Width / 2d, rect.Y + rect.Height / 2d); var pos = projection.ViewportTransform.Transform(center); - if (projection.IsCylindrical && + if (projection.IsNormalCylindrical && (pos.X < 0d || pos.X > parentMap.RenderSize.Width || pos.Y < 0d || pos.Y > parentMap.RenderSize.Height)) { diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index 41ddcf19..294cb596 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -23,6 +23,9 @@ namespace MapControl public const double PixelPerDegree = TileSize / 360d; public const double Wgs84EquatorialRadius = 6378137d; + public const double Wgs84Flattening = 1d / 298.257223563; + public static readonly double Wgs84Eccentricity = Math.Sqrt((2d - Wgs84Flattening) * Wgs84Flattening); + public const double MetersPerDegree = Wgs84EquatorialRadius * Math.PI / 180d; /// @@ -33,7 +36,7 @@ namespace MapControl /// /// Indicates if this is a normal cylindrical projection. /// - public bool IsCylindrical { get; protected set; } = false; + public bool IsNormalCylindrical { get; protected set; } = false; /// /// Indicates if this is a web mercator projection, i.e. compatible with MapTileLayer. @@ -57,12 +60,13 @@ namespace MapControl public Matrix ViewportTransform { get; private set; } /// - /// Gets the transform matrix from viewport coordinates to cartesian map coordinates (pixels). + /// Gets the transform matrix from viewport coordinates (pixels) to cartesian map coordinates. /// public Matrix InverseViewportTransform { get; private set; } /// - /// Gets the scaling factor from cartesian map coordinates to viewport coordinates. + /// Gets the scaling factor from cartesian map coordinates to viewport coordinates (pixels) + /// at the projection's point of true scale. /// public double ViewportScale { get; private set; } @@ -151,7 +155,7 @@ namespace MapControl /// public virtual string WmsQueryParameters(BoundingBox boundingBox) { - if (string.IsNullOrEmpty(CrsId) || !boundingBox.HasValidBounds) + if (string.IsNullOrEmpty(CrsId)) { return null; } diff --git a/MapControl/Shared/MapShape.cs b/MapControl/Shared/MapShape.cs index a992f776..db522ba0 100644 --- a/MapControl/Shared/MapShape.cs +++ b/MapControl/Shared/MapShape.cs @@ -106,7 +106,7 @@ namespace MapControl { var longitudeOffset = 0d; - if (parentMap.MapProjection.IsCylindrical && Location != null) + if (parentMap.MapProjection.IsNormalCylindrical && Location != null) { var viewportPosition = LocationToViewportPoint(Location); diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs index 04c52aa8..50a0c0c9 100644 --- a/MapControl/Shared/OrthographicProjection.cs +++ b/MapControl/Shared/OrthographicProjection.cs @@ -10,7 +10,7 @@ using System.Windows; namespace MapControl { /// - /// Transforms map coordinates according to the Orthographic Projection. + /// Spherical Orthographic Projection. /// public class OrthographicProjection : AzimuthalProjection { @@ -34,10 +34,11 @@ namespace MapControl var lat0 = ProjectionCenter.Latitude * Math.PI / 180d; var lat = location.Latitude * Math.PI / 180d; var dLon = (location.Longitude - ProjectionCenter.Longitude) * Math.PI / 180d; + var s = TrueScale * 180d / Math.PI; return new Point( - Wgs84EquatorialRadius * Math.Cos(lat) * Math.Sin(dLon), - Wgs84EquatorialRadius * (Math.Cos(lat0) * Math.Sin(lat) - Math.Sin(lat0) * Math.Cos(lat) * Math.Cos(dLon))); + s * Math.Cos(lat) * Math.Sin(dLon), + s * (Math.Cos(lat0) * Math.Sin(lat) - Math.Sin(lat0) * Math.Cos(lat) * Math.Cos(dLon))); } public override Location PointToLocation(Point point) @@ -47,8 +48,9 @@ namespace MapControl return ProjectionCenter; } - var x = point.X / Wgs84EquatorialRadius; - var y = point.Y / Wgs84EquatorialRadius; + var s = TrueScale * 180d / Math.PI; + var x = point.X / s; + var y = point.Y / s; var r2 = x * x + y * y; if (r2 > 1d) diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs index e77ed078..611ea7ef 100644 --- a/MapControl/Shared/StereographicProjection.cs +++ b/MapControl/Shared/StereographicProjection.cs @@ -10,7 +10,7 @@ using System.Windows; namespace MapControl { /// - /// Transforms map coordinates according to the Gnomonic Projection. + /// Spherical Stereographic Projection. /// public class StereographicProjection : AzimuthalProjection { @@ -35,7 +35,7 @@ namespace MapControl GetAzimuthDistance(ProjectionCenter, location, out azimuth, out distance); - var mapDistance = 2d * Wgs84EquatorialRadius * Math.Tan(distance / 2d); + var mapDistance = Math.Tan(distance / 2d) * TrueScale * 360d / Math.PI; return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth)); } @@ -49,7 +49,7 @@ namespace MapControl var azimuth = Math.Atan2(point.X, point.Y); var mapDistance = Math.Sqrt(point.X * point.X + point.Y * point.Y); - var distance = 2d * Math.Atan(mapDistance / (2d * Wgs84EquatorialRadius)); + var distance = 2d * Math.Atan(mapDistance / (TrueScale * 360d / Math.PI)); return GetLocation(ProjectionCenter, azimuth, distance); } diff --git a/MapControl/Shared/WebMercatorProjection.cs b/MapControl/Shared/WebMercatorProjection.cs index c97a298a..48eb7257 100644 --- a/MapControl/Shared/WebMercatorProjection.cs +++ b/MapControl/Shared/WebMercatorProjection.cs @@ -10,10 +10,8 @@ using System.Windows; namespace MapControl { /// - /// Transforms map coordinates according to the Web (or Pseudo) Mercator Projection, EPSG:3857. - /// Longitude values are transformed linearly to X values in meters, by multiplying with TrueScale. - /// Latitude values in the interval [-MaxLatitude .. MaxLatitude] are transformed to Y values in meters - /// in the interval [-R*pi .. R*pi], R=Wgs84EquatorialRadius. + /// Spherical Mercator Projection, EPSG:3857. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.41-44. /// public class WebMercatorProjection : MapProjection { @@ -25,16 +23,16 @@ namespace MapControl public WebMercatorProjection(string crsId) { CrsId = crsId; - IsCylindrical = true; + IsNormalCylindrical = true; IsWebMercator = true; MaxLatitude = YToLatitude(180d); } public override Vector GetMapScale(Location location) { - var scale = ViewportScale / Math.Cos(location.Latitude * Math.PI / 180d); + var k = 1d / Math.Cos(location.Latitude * Math.PI / 180d); // p.44 (7-3) - return new Vector(scale, scale); + return new Vector(ViewportScale * k, ViewportScale * k); } public override Point LocationToPoint(Location location) @@ -63,12 +61,12 @@ namespace MapControl return double.PositiveInfinity; } - return Math.Log(Math.Tan((latitude + 90d) * Math.PI / 360d)) / Math.PI * 180d; + return Math.Log(Math.Tan((latitude + 90d) * Math.PI / 360d)) * 180d / Math.PI; } public static double YToLatitude(double y) { - return 90d - Math.Atan(Math.Exp(-y * Math.PI / 180d)) / Math.PI * 360d; + return 90d - Math.Atan(Math.Exp(-y * Math.PI / 180d)) * 360d / Math.PI; } } } diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs index 673e8c5d..90474928 100644 --- a/MapControl/Shared/WorldMercatorProjection.cs +++ b/MapControl/Shared/WorldMercatorProjection.cs @@ -10,17 +10,12 @@ using System.Windows; namespace MapControl { /// - /// Transforms map coordinates according to the "World Mercator" Projection, EPSG:3395. - /// Longitude values are transformed linearly to X values in meters, by multiplying with TrueScale. - /// Latitude values are transformed according to the elliptical versions of the Mercator equations, - /// as shown in "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.44. + /// Elliptical Mercator Projection, EPSG:3395. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.44-45. /// public class WorldMercatorProjection : MapProjection { - public const double Wgs84Flattening = 1d / 298.257223563; - public static readonly double Wgs84Eccentricity = Math.Sqrt((2d - Wgs84Flattening) * Wgs84Flattening); - - public static double MinLatitudeDelta = 1d / Wgs84EquatorialRadius; // corresponds to 1 meter + public static double ConvergenceTolerance = 1e-6; public static int MaxIterations = 10; public WorldMercatorProjection() @@ -31,7 +26,7 @@ namespace MapControl public WorldMercatorProjection(string crsId) { CrsId = crsId; - IsCylindrical = true; + IsNormalCylindrical = true; MaxLatitude = YToLatitude(180d); } @@ -39,9 +34,9 @@ namespace MapControl { var lat = location.Latitude * Math.PI / 180d; var eSinLat = Wgs84Eccentricity * Math.Sin(lat); - var scale = ViewportScale * Math.Sqrt(1d - eSinLat * eSinLat) / Math.Cos(lat); + var k = Math.Sqrt(1d - eSinLat * eSinLat) / Math.Cos(lat); // p.44 (7-8) - return new Vector(scale, scale); + return new Vector(ViewportScale * k, ViewportScale * k); } public override Point LocationToPoint(Location location) @@ -72,24 +67,24 @@ namespace MapControl var lat = latitude * Math.PI / 180d; - return Math.Log(Math.Tan(lat / 2d + Math.PI / 4d) * ConformalFactor(lat)) / Math.PI * 180d; + return Math.Log(Math.Tan(lat / 2d + Math.PI / 4d) + * ConformalFactor(lat)) * 180d / Math.PI; // p.44 (7-7) } public static double YToLatitude(double y) { - var t = Math.Exp(-y * Math.PI / 180d); - var lat = Math.PI / 2d - 2d * Math.Atan(t); - var latDelta = 1d; + var t = Math.Exp(-y * Math.PI / 180d); // p.44 (7-10) + var lat = Math.PI / 2d - 2d * Math.Atan(t); // p.44 (7-11) + var relChange = 1d; - for (int i = 0; i < MaxIterations && latDelta > MinLatitudeDelta; i++) + for (var i = 0; i < MaxIterations && relChange > ConvergenceTolerance; i++) { - var newLat = Math.PI / 2d - 2d * Math.Atan(t * ConformalFactor(lat)); - - latDelta = Math.Abs(newLat - lat); + var newLat = Math.PI / 2d - 2d * Math.Atan(t * ConformalFactor(lat)); // p.44 (7-9) + relChange = Math.Abs(1d - newLat / lat); lat = newLat; } - return lat / Math.PI * 180d; + return lat * 180d / Math.PI; } private static double ConformalFactor(double lat) diff --git a/MapControl/UWP/MapGraticule.UWP.cs b/MapControl/UWP/MapGraticule.UWP.cs index f9de665e..d939718f 100644 --- a/MapControl/UWP/MapGraticule.UWP.cs +++ b/MapControl/UWP/MapGraticule.UWP.cs @@ -26,7 +26,7 @@ namespace MapControl { var projection = ParentMap.MapProjection; - if (projection.IsCylindrical) + if (projection.IsNormalCylindrical) { if (path == null) { diff --git a/MapControl/WPF/MapGraticule.WPF.cs b/MapControl/WPF/MapGraticule.WPF.cs index 58727bfe..4b16c773 100644 --- a/MapControl/WPF/MapGraticule.WPF.cs +++ b/MapControl/WPF/MapGraticule.WPF.cs @@ -44,7 +44,7 @@ namespace MapControl var lineDistance = GetLineDistance(); var labelFormat = GetLabelFormat(lineDistance); - if (projection.IsCylindrical) + if (projection.IsNormalCylindrical) { DrawCylindricalGraticule(drawingContext, projection, lineDistance, labelFormat); } @@ -58,52 +58,48 @@ namespace MapControl private void DrawCylindricalGraticule(DrawingContext drawingContext, MapProjection projection, double lineDistance, string labelFormat) { var boundingBox = projection.ViewportRectToBoundingBox(new Rect(ParentMap.RenderSize)); + var latLabelStart = Math.Ceiling(boundingBox.South / lineDistance) * lineDistance; + var lonLabelStart = Math.Ceiling(boundingBox.West / lineDistance) * lineDistance; + var latLabels = new List public class GeoApiProjection : MapProjection { - private ICoordinateTransformation coordinateTransform; - private IMathTransform mathTransform; - private IMathTransform inverseTransform; + private IProjectedCoordinateSystem coordinateSystem; + + public IMathTransform MathTransform { get; private set; } + public IMathTransform InverseTransform { get; private set; } /// - /// Gets or sets the underlying ICoordinateTransformation instance. - /// Setting this property updates the CrsId property. + /// Gets or sets the IProjectedCoordinateSystem of the MapProjection. /// - public ICoordinateTransformation CoordinateTransform + public IProjectedCoordinateSystem CoordinateSystem { - get { return coordinateTransform; } + get { return coordinateSystem; } set { if (value == null) @@ -39,15 +37,45 @@ namespace MapControl.Projections throw new ArgumentNullException("The property value must not be null."); } - coordinateTransform = value; - mathTransform = coordinateTransform.MathTransform; - inverseTransform = mathTransform.Inverse(); + coordinateSystem = value; - if (coordinateTransform.TargetCS != null && - !string.IsNullOrEmpty(coordinateTransform.TargetCS.Authority) && - coordinateTransform.TargetCS.AuthorityCode > 0) + var coordinateTransform = new CoordinateTransformationFactory() + .CreateFromCoordinateSystems(GeographicCoordinateSystem.WGS84, coordinateSystem); + + MathTransform = coordinateTransform.MathTransform; + InverseTransform = MathTransform.Inverse(); + + CrsId = (!string.IsNullOrEmpty(coordinateSystem.Authority) && coordinateSystem.AuthorityCode > 0) + ? string.Format("{0}:{1}", coordinateSystem.Authority, coordinateSystem.AuthorityCode) + : null; + + if (!IsWebMercator) { - CrsId = string.Format("{0}:{1}", coordinateTransform.TargetCS.Authority, coordinateTransform.TargetCS.AuthorityCode); + IsWebMercator = CrsId == "EPSG:3857" || CrsId == "EPSG:900913"; + } + + var projection = coordinateSystem.Projection; + var scaleFactor = projection.GetParameter("scale_factor"); + + if (scaleFactor != null) + { + TrueScale = scaleFactor.Value * MetersPerDegree; + } + + if (!IsNormalCylindrical) + { + var centralMeridian = projection.GetParameter("central_meridian") ?? projection.GetParameter("longitude_of_origin"); + var centralParallel = projection.GetParameter("latitude_of_origin") ?? projection.GetParameter("central_parallel"); + var falseEasting = projection.GetParameter("false_easting"); + var falseNorthing = projection.GetParameter("false_northing"); + + if (centralMeridian != null && centralMeridian.Value == 0d && + centralParallel != null && centralParallel.Value == 0d && + (falseEasting == null || falseEasting.Value == 0d) && + (falseNorthing == null || falseNorthing.Value == 0d)) + { + IsNormalCylindrical = true; + } } } } @@ -55,40 +83,34 @@ namespace MapControl.Projections /// /// Gets or sets an OGC Well-known text representation of a projected coordinate system, /// i.e. a PROJCS[...] string as used by https://epsg.io or http://spatialreference.org. - /// Setting this property updates the CoordinateTransform property. + /// Setting this property updates the CoordinateSystem property with an IProjectedCoordinateSystem created from the WKT string. /// public string WKT { - get { return coordinateTransform?.TargetCS?.WKT; } - set - { - var sourceCs = GeographicCoordinateSystem.WGS84; - var targetCs = (ICoordinateSystem)CoordinateSystemWktReader.Parse(value, Encoding.UTF8); - - CoordinateTransform = new CoordinateTransformationFactory().CreateFromCoordinateSystems(sourceCs, targetCs); - } + get { return CoordinateSystem?.WKT; } + set { CoordinateSystem = (IProjectedCoordinateSystem)new CoordinateSystemFactory().CreateFromWkt(value); } } public override Point LocationToPoint(Location location) { - if (mathTransform == null) + if (MathTransform == null) { - throw new InvalidOperationException("The CoordinateTransformation property is not set."); + throw new InvalidOperationException("The CoordinateSystem property is not set."); } - var coordinate = mathTransform.Transform(new Coordinate(location.Longitude, location.Latitude)); + var coordinate = MathTransform.Transform(new Coordinate(location.Longitude, location.Latitude)); return new Point(coordinate.X, coordinate.Y); } public override Location PointToLocation(Point point) { - if (inverseTransform == null) + if (InverseTransform == null) { - throw new InvalidOperationException("The CoordinateTransformation property is not set."); + throw new InvalidOperationException("The CoordinateSystem property is not set."); } - var coordinate = inverseTransform.Transform(new Coordinate(point.X, point.Y)); + var coordinate = InverseTransform.Transform(new Coordinate(point.X, point.Y)); return new Location(coordinate.Y, coordinate.X); } diff --git a/MapProjections/Shared/PolarStereographicProjection.cs b/MapProjections/Shared/PolarStereographicProjection.cs new file mode 100644 index 00000000..a91ff867 --- /dev/null +++ b/MapProjections/Shared/PolarStereographicProjection.cs @@ -0,0 +1,129 @@ +// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control +// © 2018 Clemens Fischer +// Licensed under the Microsoft Public License (Ms-PL) + +using System; +#if !WINDOWS_UWP +using System.Windows; +#endif + +namespace MapControl.Projections +{ + /// + /// Elliptical Polar Stereographic Projection with a given scale factor at the pole and + /// optional false easting and northing, as used by the UPS North and UPS South projections. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.160-162. + /// + public class PolarStereographicProjection : MapProjection + { + public static double ConvergenceTolerance = 1e-6; + public static int MaxIterations = 10; + + private readonly bool north; + private readonly double scaleFactor; + private readonly double falseEasting; + private readonly double falseNorthing; + + public PolarStereographicProjection(string crsId, bool north, double scaleFactor = 1d, double falseEasting = 0d, double falseNorthing = 0d) + { + CrsId = crsId; + TrueScale = scaleFactor * MetersPerDegree; + this.north = north; + this.scaleFactor = scaleFactor; + this.falseEasting = falseEasting; + this.falseNorthing = falseNorthing; + } + + public override Vector GetMapScale(Location location) + { + var lat = (north ? location.Latitude : -location.Latitude) * Math.PI / 180d; + var a = Wgs84EquatorialRadius; + var e = Wgs84Eccentricity; + var s = Math.Sqrt(Math.Pow(1 + e, 1 + e) * Math.Pow(1 - e, 1 - e)); + var t = Math.Tan(Math.PI / 4d - lat / 2d) / ConformalFactor(lat); + var rho = 2d * a * scaleFactor * t / s; + var eSinLat = e * Math.Sin(lat); + var m = Math.Cos(lat) / Math.Sqrt(1d - eSinLat * eSinLat); + var k = rho / (a * m); + + return new Vector(ViewportScale * k, ViewportScale * k); + } + + public override Point LocationToPoint(Location location) + { + var lat = location.Latitude * Math.PI / 180d; + var lon = location.Longitude * Math.PI / 180d; + + if (north) + { + lon = Math.PI - lon; + } + else + { + lat = -lat; + } + + var a = Wgs84EquatorialRadius; + var e = Wgs84Eccentricity; + var s = Math.Sqrt(Math.Pow(1 + e, 1 + e) * Math.Pow(1 - e, 1 - e)); + var t = Math.Tan(Math.PI / 4d - lat / 2d) / ConformalFactor(lat); + var rho = 2d * a * scaleFactor * t / s; + + return new Point(rho * Math.Sin(lon) + falseEasting, rho * Math.Cos(lon) + falseNorthing); + } + + public override Location PointToLocation(Point point) + { + point.X -= falseEasting; + point.Y -= falseNorthing; + + var lon = Math.Atan2(point.X, point.Y); + var rho = Math.Sqrt(point.X * point.X + point.Y * point.Y); + var a = Wgs84EquatorialRadius; + var e = Wgs84Eccentricity; + var s = Math.Sqrt(Math.Pow(1 + e, 1 + e) * Math.Pow(1 - e, 1 - e)); + var t = rho * s / (2d * a * scaleFactor); + var lat = Math.PI / 2d - 2d * Math.Atan(t); + var relChange = 1d; + + for (int i = 0; i < MaxIterations && relChange > ConvergenceTolerance; i++) + { + var newLat = Math.PI / 2d - 2d * Math.Atan(t * ConformalFactor(lat)); + relChange = Math.Abs(1d - newLat / lat); + lat = newLat; + } + + if (north) + { + lon = Math.PI - lon; + } + else + { + lat = -lat; + } + + return new Location(lat * 180d / Math.PI, lon * 180d / Math.PI); + } + + private static double ConformalFactor(double lat) + { + var eSinLat = Wgs84Eccentricity * Math.Sin(lat); + + return Math.Pow((1d - eSinLat) / (1d + eSinLat), Wgs84Eccentricity / 2d); + } + } + + public class UpsNorthProjection : PolarStereographicProjection + { + public UpsNorthProjection() : base("EPSG:32661", true, 0.994, 2e6, 2e6) + { + } + } + + public class UpsSouthProjection : PolarStereographicProjection + { + public UpsSouthProjection() : base("EPSG:32761", false, 0.994, 2e6, 2e6) + { + } + } +} diff --git a/MapProjections/Shared/UtmProjection.cs b/MapProjections/Shared/UtmProjection.cs index d1856e9c..ccd2a7c3 100644 --- a/MapProjections/Shared/UtmProjection.cs +++ b/MapProjections/Shared/UtmProjection.cs @@ -3,21 +3,12 @@ // Licensed under the Microsoft Public License (Ms-PL) using System; -#if !WINDOWS_UWP -using System.Windows; -#endif +using ProjNet.CoordinateSystems; namespace MapControl.Projections { public class UtmProjection : GeoApiProjection { - private const string wktFormat = "PROJCS[\"WGS 84 / UTM zone {0}\", GEOGCS[\"WGS 84\", DATUM[\"WGS_1984\", SPHEROID[\"WGS 84\", 6378137, 298.257223563, AUTHORITY[\"EPSG\", \"7030\"]], AUTHORITY[\"EPSG\", \"6326\"]], PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]], UNIT[\"degree\", 0.01745329251994328, AUTHORITY[\"EPSG\", \"9122\"]], AUTHORITY[\"EPSG\", \"4326\"]], UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]], PROJECTION[\"Transverse_Mercator\"], PARAMETER[\"latitude_of_origin\", 0], PARAMETER[\"central_meridian\", {1}], PARAMETER[\"scale_factor\", 0.9996], PARAMETER[\"false_easting\", 500000], PARAMETER[\"false_northing\", {2}], AUTHORITY[\"EPSG\", \"{3}\"], AXIS[\"Easting\", EAST], AXIS[\"Northing\", NORTH]]"; - - public UtmProjection() - { - TrueScale = 0.9996 * MetersPerDegree; - } - private string zone; public string Zone @@ -48,7 +39,7 @@ namespace MapControl.Projections public void SetZone(int zoneNumber, bool north) { - if (zoneNumber < 1 || zoneNumber > 60) + if (zoneNumber < 1 || zoneNumber > 61) { throw new ArgumentException("Invalid UTM zone number."); } @@ -57,12 +48,8 @@ namespace MapControl.Projections if (zone != zoneName) { - var centralMeridian = zoneNumber * 6 - 183; - var falseNorthing = north ? 0 : 10000000; - var authorityCode = (north ? 32600 : 32700) + zoneNumber; - zone = zoneName; - WKT = string.Format(wktFormat, zone, centralMeridian, falseNorthing, authorityCode); + CoordinateSystem = ProjectedCoordinateSystem.WGS84_UTM(zoneNumber, north); } } diff --git a/MapProjections/Shared/WebMercatorProjection.cs b/MapProjections/Shared/WebMercatorProjection.cs new file mode 100644 index 00000000..314eef3d --- /dev/null +++ b/MapProjections/Shared/WebMercatorProjection.cs @@ -0,0 +1,33 @@ +// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control +// © 2018 Clemens Fischer +// Licensed under the Microsoft Public License (Ms-PL) + +using System; +#if !WINDOWS_UWP +using System.Windows; +#endif +using ProjNet.CoordinateSystems; + +namespace MapControl.Projections +{ + /// + /// Spherical Mercator Projection implemented by setting the CoordinateSystem property of a GeoApiProjection. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.41-44. + /// + public class WebMercatorProjection : GeoApiProjection + { + public WebMercatorProjection() + { + IsWebMercator = true; + IsNormalCylindrical = true; + CoordinateSystem = ProjectedCoordinateSystem.WebMercator; + } + + public override Vector GetMapScale(Location location) + { + var k = 1d / Math.Cos(location.Latitude * Math.PI / 180d); // p.44 (7-3) + + return new Vector(ViewportScale * k, ViewportScale * k); + } + } +} diff --git a/MapProjections/Shared/WorldMercatorProjection.cs b/MapProjections/Shared/WorldMercatorProjection.cs new file mode 100644 index 00000000..ada2428c --- /dev/null +++ b/MapProjections/Shared/WorldMercatorProjection.cs @@ -0,0 +1,48 @@ +// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control +// © 2018 Clemens Fischer +// Licensed under the Microsoft Public License (Ms-PL) + +using System; +#if !WINDOWS_UWP +using System.Windows; +#endif + +namespace MapControl.Projections +{ + /// + /// Elliptical Mercator Projection implemented by setting the WKT property of a GeoApiProjection. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.44-45. + /// + public class WorldMercatorProjection : GeoApiProjection + { + public WorldMercatorProjection() + { + IsNormalCylindrical = true; + WKT = "PROJCS[\"WGS 84 / World Mercator\"," + + "GEOGCS[\"WGS 84\"," + + "DATUM[\"WGS_1984\", SPHEROID[\"WGS 84\", 6378137, 298.257223563, AUTHORITY[\"EPSG\", \"7030\"]], AUTHORITY[\"EPSG\", \"6326\"]]," + + "PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]]," + + "UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9122\"]]," + + "AUTHORITY[\"EPSG\", \"4326\"]]," + + "PROJECTION[\"Mercator_1SP\"]," + + "PARAMETER[\"latitude_of_origin\", 0]," + + "PARAMETER[\"central_meridian\", 0]," + + "PARAMETER[\"scale_factor\", 1]," + + "PARAMETER[\"false_easting\", 0]," + + "PARAMETER[\"false_northing\", 0]," + + "UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]]," + + "AXIS[\"Easting\", EAST]," + + "AXIS[\"Northing\", NORTH]," + + "AUTHORITY[\"EPSG\", \"3395\"]]"; + } + + public override Vector GetMapScale(Location location) + { + var lat = location.Latitude * Math.PI / 180d; + var eSinLat = Wgs84Eccentricity * Math.Sin(lat); + var k = Math.Sqrt(1d - eSinLat * eSinLat) / Math.Cos(lat); // p.44 (7-8) + + return new Vector(ViewportScale * k, ViewportScale * k); + } + } +} diff --git a/MapProjections/UWP/MapProjections.UWP.csproj b/MapProjections/UWP/MapProjections.UWP.csproj index 6424f1c2..dbe705a6 100644 --- a/MapProjections/UWP/MapProjections.UWP.csproj +++ b/MapProjections/UWP/MapProjections.UWP.csproj @@ -43,9 +43,15 @@ GeoApiProjection.cs + + PolarStereographicProjection.cs + UtmProjection.cs + + WebMercatorProjection.cs + diff --git a/MapProjections/WPF/MapProjections.WPF.csproj b/MapProjections/WPF/MapProjections.WPF.csproj index ccaa6793..5d92952f 100644 --- a/MapProjections/WPF/MapProjections.WPF.csproj +++ b/MapProjections/WPF/MapProjections.WPF.csproj @@ -56,9 +56,18 @@ GeoApiProjection.cs + + PolarStereographicProjection.cs + UtmProjection.cs + + WebMercatorProjection.cs + + + WorldMercatorProjection.cs + Code