Updated PolarStereographicProjection

This commit is contained in:
ClemensFischer 2026-01-18 21:03:25 +01:00
parent 31b38d5049
commit 19cba8978b
3 changed files with 13 additions and 19 deletions

View file

@ -112,24 +112,16 @@ namespace MapControl
{
var lat1 = Latitude * Math.PI / 180d;
var lon1 = Longitude * Math.PI / 180d;
var cosD = Math.Cos(distance);
var sinD = Math.Sin(distance);
var cosA = Math.Cos(azimuth);
var sinA = Math.Sin(azimuth);
var cosLat1 = Math.Cos(lat1);
var sinLat1 = Math.Sin(lat1);
var cosA = Math.Cos(azimuth);
var sinA = Math.Sin(azimuth);
var cosD = Math.Cos(distance);
var sinD = Math.Sin(distance);
var lat2 = Math.Asin(sinLat1 * cosD + cosLat1 * sinD * cosA);
var lon2 = lon1 + Math.Atan2(sinD * sinA, cosLat1 * cosD - sinLat1 * sinD * cosA);
return new Location(lat2 * 180d / Math.PI, lon2 * 180d / Math.PI);
}
/// <summary>
/// Calculates the Location on a great circle at the specified azimuth in degrees and distance in meters from this Location.
/// </summary>
public Location GetLocation(double azimuth, double distance, double earthRadius = MapProjection.Wgs84MeanRadius)
{
return GetLocation(azimuth * Math.PI / 180d, distance / earthRadius);
}
}
}

View file

@ -26,7 +26,7 @@ namespace MapControl
public double FalseNorthing { get; set; } = 2e6;
public Hemisphere Hemisphere { get; set; }
public static double RelativeScale(Hemisphere hemisphere, double equatorialRadius, double flattening, double scaleFactor, double latitude)
public static double RelativeScale(Hemisphere hemisphere, double flattening, double scaleFactor, double latitude)
{
var sign = hemisphere == Hemisphere.North ? 1d : -1d;
var phi = sign * latitude * Math.PI / 180d;
@ -36,17 +36,19 @@ namespace MapControl
var t = Math.Tan(Math.PI / 4d - phi / 2d)
/ Math.Pow((1d - eSinPhi) / (1d + eSinPhi), e / 2d); // p.161 (15-9)
var r = 2d * equatorialRadius * scaleFactor * t
// r == ρ/a
var r = 2d * scaleFactor * t
/ Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e)); // p.161 (21-33)
var m = Math.Cos(phi) / Math.Sqrt(1d - eSinPhi * eSinPhi); // p.160 (14-15)
return r / (equatorialRadius * m); // p.161 (21-32)
return r / m; // p.161 (21-32)
}
public override Point RelativeScale(double latitude, double longitude)
{
var k = RelativeScale(Hemisphere, EquatorialRadius, Flattening, ScaleFactor, latitude);
var k = RelativeScale(Hemisphere, Flattening, ScaleFactor, latitude);
return new Point(k, k);
}
@ -62,7 +64,7 @@ namespace MapControl
var t = Math.Tan(Math.PI / 4d - phi / 2d)
/ Math.Pow((1d - eSinPhi) / (1d + eSinPhi), e / 2d); // p.161 (15-9)
// ρ
var r = 2d * EquatorialRadius * ScaleFactor * t
/ Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e)); // p.161 (21-33)

View file

@ -25,7 +25,7 @@ namespace MapControl.Projections
public override Point RelativeScale(double latitude, double longitude)
{
var k = PolarStereographicProjection.RelativeScale(Hemisphere.North, Wgs84EquatorialRadius, Wgs84Flattening, 0.994, latitude);
var k = PolarStereographicProjection.RelativeScale(Hemisphere.North, Wgs84Flattening, 0.994, latitude);
return new Point(k, k);
}
@ -50,7 +50,7 @@ namespace MapControl.Projections
public override Point RelativeScale(double latitude, double longitude)
{
var k = PolarStereographicProjection.RelativeScale(Hemisphere.South, Wgs84EquatorialRadius, Wgs84Flattening, 0.994, latitude);
var k = PolarStereographicProjection.RelativeScale(Hemisphere.South, Wgs84Flattening, 0.994, latitude);
return new Point(k, k);
}