Updated projections

This commit is contained in:
ClemensFischer 2026-01-13 23:24:48 +01:00
parent da3a9a18fb
commit 91dfeee19c
6 changed files with 49 additions and 43 deletions

View file

@ -32,18 +32,19 @@ namespace MapControl
return new Point();
}
var lat0 = Center.Latitude * Math.PI / 180d;
var lat = latitude * Math.PI / 180d;
var dLon = (longitude - Center.Longitude) * Math.PI / 180d;
var phi = latitude * Math.PI / 180d;
var phi1 = Center.Latitude * Math.PI / 180d;
var lambda = (longitude - Center.Longitude) * Math.PI / 180d; // λ - λ0
if (Math.Abs(lat - lat0) > Math.PI / 2d || Math.Abs(dLon) > Math.PI / 2d)
if (Math.Abs(phi - phi1) > Math.PI / 2d || Math.Abs(lambda) > Math.PI / 2d)
{
return null;
}
return new Point(
Wgs84MeanRadius * Math.Cos(lat) * Math.Sin(dLon),
Wgs84MeanRadius * (Math.Cos(lat0) * Math.Sin(lat) - Math.Sin(lat0) * Math.Cos(lat) * Math.Cos(dLon)));
var x = Wgs84MeanRadius * Math.Cos(phi) * Math.Sin(lambda); // p.149 (20-3)
var y = Wgs84MeanRadius * (Math.Cos(phi1) * Math.Sin(phi) -
Math.Sin(phi1) * Math.Cos(phi) * Math.Cos(lambda)); // p.149 (20-4)
return new Point(x, y);
}
public override Location MapToLocation(double x, double y)
@ -66,13 +67,14 @@ namespace MapControl
var sinC = r;
var cosC = Math.Sqrt(1 - r2);
var lat0 = Center.Latitude * Math.PI / 180d;
var cosLat0 = Math.Cos(lat0);
var sinLat0 = Math.Sin(lat0);
var phi1 = Center.Latitude * Math.PI / 180d;
var cosPhi1 = Math.Cos(phi1);
var sinPhi1 = Math.Sin(phi1);
return new Location(
180d / Math.PI * Math.Asin(cosC * sinLat0 + y * sinC * cosLat0 / r),
180d / Math.PI * Math.Atan2(x * sinC, r * cosC * cosLat0 - y * sinC * sinLat0) + Center.Longitude);
var phi = Math.Asin(cosC * sinPhi1 + y * sinC * cosPhi1 / r); // p.150 (20-14)
var lambda = Math.Atan2(x * sinC, r * cosC * cosPhi1 - y * sinC * sinPhi1); // p.150 (20-15)
return new Location(180d / Math.PI * phi, 180d / Math.PI * lambda + Center.Longitude);
}
}
}

View file

@ -96,16 +96,16 @@ namespace MapControl
var t = r * Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e))
/ (2d * EquatorialRadius * ScaleFactor); // p.162 (21-39)
var lat = WorldMercatorProjection.ApproximateLatitude(e, t); // p.162 (3-5)
var lon = Math.Atan2(x, -y); // p.162 (20-16)
var phi = WorldMercatorProjection.ApproximateLatitude(e, t); // p.162 (3-5)
var lambda = Math.Atan2(x, -y); // p.162 (20-16)
if (Hemisphere == Hemisphere.South)
{
lat = -lat;
lon = -lon;
phi = -phi;
lambda = -lambda;
}
return new Location(lat * 180d / Math.PI, lon * 180d / Math.PI);
return new Location(phi * 180d / Math.PI, lambda * 180d / Math.PI);
}
}

View file

@ -30,9 +30,7 @@ namespace MapControl
public override Point RelativeScale(double latitude, double longitude)
{
var lat = latitude * Math.PI / 180d;
var eSinLat = Wgs84Eccentricity * Math.Sin(lat);
var k = Math.Sqrt(1d - eSinLat * eSinLat) / Math.Cos(lat); // p.44 (7-8)
var k = ScaleFactor(latitude);
return new Point(k, k);
}
@ -51,6 +49,14 @@ namespace MapControl
x / Wgs84MeterPerDegree);
}
public static double ScaleFactor(double latitude)
{
var phi = latitude * Math.PI / 180d;
var eSinPhi = Wgs84Eccentricity * Math.Sin(phi);
return Math.Sqrt(1d - eSinPhi * eSinPhi) / Math.Cos(phi); // p.44 (7-8)
}
public static double LatitudeToY(double latitude)
{
if (latitude <= -90d)
@ -63,11 +69,11 @@ namespace MapControl
return double.PositiveInfinity;
}
var lat = latitude * Math.PI / 180d;
var eSinLat = Wgs84Eccentricity * Math.Sin(lat);
var f = Math.Pow((1d - eSinLat) / (1d + eSinLat), Wgs84Eccentricity / 2d);
var phi = latitude * Math.PI / 180d;
var eSinPhi = Wgs84Eccentricity * Math.Sin(phi);
var f = Math.Pow((1d - eSinPhi) / (1d + eSinPhi), Wgs84Eccentricity / 2d);
return Math.Log(Math.Tan(lat / 2d + Math.PI / 4d) * f) * 180d / Math.PI; // p.44 (7-7)
return Math.Log(Math.Tan(phi / 2d + Math.PI / 4d) * f) * 180d / Math.PI; // p.44 (7-7)
}
public static double YToLatitude(double y)
@ -79,18 +85,17 @@ namespace MapControl
internal static double ApproximateLatitude(double e, double t)
{
var e_2 = e * e;
var e_4 = e_2 * e_2;
var e_6 = e_2 * e_4;
var e_8 = e_2 * e_6;
var e2 = e * e;
var e4 = e2 * e2;
var e6 = e2 * e4;
var e8 = e2 * e6;
var chi = Math.PI / 2d - 2d * Math.Atan(t); // p.45 (7-13)
var lat = Math.PI / 2d - 2d * Math.Atan(t); // p.45 (7-13)
return lat
+ (e_2 / 2d + 5d * e_4 / 24d + e_6 / 12d + 13d * e_8 / 360d) * Math.Sin(2d * lat)
+ (7d * e_4 / 48d + 29d * e_6 / 240d + 811d * e_8 / 11520d) * Math.Sin(4d * lat)
+ (7d * e_6 / 120d + 81d * e_8 / 1120d) * Math.Sin(6d * lat)
+ 4279d * e_8 / 161280d * Math.Sin(8d * lat); // p.45 (3-5)
return chi
+ (e2 / 2d + 5d * e4 / 24d + e6 / 12d + 13d * e8 / 360d) * Math.Sin(2d * chi)
+ (7d * e4 / 48d + 29d * e6 / 240d + 811d * e8 / 11520d) * Math.Sin(4d * chi)
+ (7d * e6 / 120d + 81d * e8 / 1120d) * Math.Sin(6d * chi)
+ 4279d * e8 / 161280d * Math.Sin(8d * chi); // p.45 (3-5)
}
}
}

View file

@ -241,7 +241,7 @@ namespace MapControl.Projections
{
MapControl.WebMercatorProjection.DefaultCrsId => new WebMercatorProjection(),
MapControl.WorldMercatorProjection.DefaultCrsId => new WorldMercatorProjection(),
MapControl.Wgs84AutoUtmProjection.DefaultCrsId => new Wgs84AutoUtmProjection(),
Wgs84AutoUtmProjection.DefaultCrsId => new Wgs84AutoUtmProjection(),
_ => base.GetProjection(crsId)
};

View file

@ -9,10 +9,12 @@ namespace MapControl.Projections
/// </summary>
public class Wgs84AutoUtmProjection : Wgs84UtmProjection
{
public const string DefaultCrsId = "AUTO2:42001";
private readonly string autoCrsId;
public Wgs84AutoUtmProjection() // parameterless constructor for XAML
: this(MapControl.Wgs84AutoUtmProjection.DefaultCrsId)
: this(DefaultCrsId)
{
}

View file

@ -1,5 +1,4 @@
using System;
#if WPF
#if WPF
using System.Windows;
#elif AVALONIA
using Avalonia;
@ -36,9 +35,7 @@ namespace MapControl.Projections
public override Point RelativeScale(double latitude, double longitude)
{
var lat = latitude * Math.PI / 180d;
var eSinLat = MapControl.WorldMercatorProjection.Wgs84Eccentricity * Math.Sin(lat);
var k = Math.Sqrt(1d - eSinLat * eSinLat) / Math.Cos(lat); // p.44 (7-8)
var k = MapControl.WorldMercatorProjection.ScaleFactor(latitude);
return new Point(k, k);
}