Rotated relative scale in MapProjections

This commit is contained in:
ClemensFischer 2026-01-20 15:48:30 +01:00
parent 25d4e13c16
commit 3fbfb0d5c1
8 changed files with 85 additions and 41 deletions

View file

@ -28,12 +28,12 @@ namespace MapControl
public override Matrix RelativeScale(double latitude, double longitude)
{
(var cosC, var _, var _) = GetPointValues(latitude, longitude);
var p = GetProjectedPoint(latitude, longitude);
var k = 1d;
if (cosC < 1d)
if (p.CosC < 1d)
{
var c = Math.Acos(cosC);
var c = Math.Acos(p.CosC);
k = c / Math.Sin(c); // p.195 (25-2)
}
@ -42,16 +42,16 @@ namespace MapControl
public override Point? LocationToMap(double latitude, double longitude)
{
(var cosC, var x, var y) = GetPointValues(latitude, longitude);
var p = GetProjectedPoint(latitude, longitude);
var k = 1d;
if (cosC < 1d)
if (p.CosC < 1d)
{
var c = Math.Acos(cosC);
var c = Math.Acos(p.CosC);
k = c / Math.Sin(c); // p.195 (25-2)
}
return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.195 (22-4/5)
return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.195 (22-4/5)
}
public override Location MapToLocation(double x, double y)

View file

@ -3,7 +3,7 @@
namespace MapControl
{
/// <summary>
/// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395).
/// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.141.
/// </summary>
public abstract class AzimuthalProjection : MapProjection
{
@ -14,22 +14,33 @@ namespace MapControl
public double EarthRadius { get; set; } = Wgs84MeanRadius;
protected (double, double, double) GetPointValues(double latitude, double longitude)
public readonly struct ProjectedPoint
{
var phi = latitude * Math.PI / 180d;
var phi1 = Center.Latitude * Math.PI / 180d;
var dLambda = (longitude - Center.Longitude) * Math.PI / 180d; // λ - λ0
var cosPhi = Math.Cos(phi);
var sinPhi = Math.Sin(phi);
var cosPhi1 = Math.Cos(phi1);
var sinPhi1 = Math.Sin(phi1);
var cosLambda = Math.Cos(dLambda);
var sinLambda = Math.Sin(dLambda);
var cosC = sinPhi1 * sinPhi + cosPhi1 * cosPhi * cosLambda; // (5-3)
var x = cosPhi * sinLambda;
var y = cosPhi1 * sinPhi - sinPhi1 * cosPhi * cosLambda;
public double X { get; }
public double Y { get; }
public double CosC { get; }
return (cosC, x, y);
public ProjectedPoint(double centerLatitude, double centerLongitude, double latitude, double longitude)
{
var phi = latitude * Math.PI / 180d;
var phi1 = centerLatitude * Math.PI / 180d;
var dLambda = (longitude - centerLongitude) * Math.PI / 180d; // λ - λ0
var cosPhi = Math.Cos(phi);
var sinPhi = Math.Sin(phi);
var cosPhi1 = Math.Cos(phi1);
var sinPhi1 = Math.Sin(phi1);
var cosLambda = Math.Cos(dLambda);
var sinLambda = Math.Sin(dLambda);
X = cosPhi * sinLambda;
Y = cosPhi1 * sinPhi - sinPhi1 * cosPhi * cosLambda;
CosC = sinPhi1 * sinPhi + cosPhi1 * cosPhi * cosLambda; // (5-3)
}
}
protected ProjectedPoint GetProjectedPoint(double latitude, double longitude)
{
return new ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude);
}
protected Location GetLocation(double x, double y, double rho, double sinC)

View file

@ -28,27 +28,27 @@ namespace MapControl
public override Matrix RelativeScale(double latitude, double longitude)
{
(var cosC, var x, var y) = GetPointValues(latitude, longitude);
var k = 1d / cosC; // p.165 (22-3)
var p = GetProjectedPoint(latitude, longitude);
var k = 1d / p.CosC; // p.165 (22-3)
var h = k * k; // p.165 (22-2)
var scale = new Matrix(h, 0d, 0d, k, 0d, 0d);
scale.Rotate(-Math.Atan2(y, x) * 180d / Math.PI);
scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI);
return scale;
}
public override Point? LocationToMap(double latitude, double longitude)
{
(var cosC, var x, var y) = GetPointValues(latitude, longitude);
var p = GetProjectedPoint(latitude, longitude);
if (cosC <= 0d)
if (p.CosC <= 0d)
{
return null;
}
var k = 1d / cosC; // p.165 (22-3)
var k = 1d / p.CosC; // p.165 (22-3)
return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.165 (22-4/5)
return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.165 (22-4/5)
}
public override Location MapToLocation(double x, double y)

View file

@ -28,8 +28,8 @@ namespace MapControl
public abstract class MapProjection
{
public const double Wgs84EquatorialRadius = 6378137d;
public const double Wgs84MeterPerDegree = Wgs84EquatorialRadius * Math.PI / 180d;
public const double Wgs84Flattening = 1d / 298.257223563;
public const double Wgs84MeterPerDegree = Wgs84EquatorialRadius * Math.PI / 180d;
// Arithmetic mean radius (2*a + b) / 3 == (1 - f/3) * a.
// See https://en.wikipedia.org/wiki/Earth_radius#Arithmetic_mean_radius.

View file

@ -28,19 +28,19 @@ namespace MapControl
public override Matrix RelativeScale(double latitude, double longitude)
{
(var cosC, var x, var y) = GetPointValues(latitude, longitude);
var h = cosC; // p.149 (20-5)
var p = GetProjectedPoint(latitude, longitude);
var h = p.CosC; // p.149 (20-5)
var scale = new Matrix(h, 0d, 0d, 1d, 0d, 0d);
scale.Rotate(-Math.Atan2(y, x) * 180d / Math.PI);
scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI);
return scale;
}
public override Point? LocationToMap(double latitude, double longitude)
{
(var cosC, var x, var y) = GetPointValues(latitude, longitude);
var p = GetProjectedPoint(latitude, longitude);
return cosC >= 0d ? new Point(EarthRadius * x, EarthRadius * y) : null; // p.149 (20-3/4)
return p.CosC >= 0d ? new Point(EarthRadius * p.X, EarthRadius * p.Y) : null; // p.149 (20-3/4)
}
public override Location MapToLocation(double x, double y)

View file

@ -28,18 +28,18 @@ namespace MapControl
public override Matrix RelativeScale(double latitude, double longitude)
{
(var cosC, var _, var _) = GetPointValues(latitude, longitude);
var k = 2d / (1d + cosC); // p.157 (21-4), k0 == 1
var p = GetProjectedPoint(latitude, longitude);
var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1
return new Matrix(k, 0d, 0d, k, 0d, 0d);
}
public override Point? LocationToMap(double latitude, double longitude)
{
(var cosC, var x, var y) = GetPointValues(latitude, longitude);
var k = 2d / (1d + cosC); // p.157 (21-4), k0 == 1
var p = GetProjectedPoint(latitude, longitude);
var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1
return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.157 (21-2/3)
return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.157 (21-2/3)
}
public override Location MapToLocation(double x, double y)

View file

@ -1,7 +1,15 @@
using System.Globalization;
using System;
using System.Globalization;
#if WPF
using System.Windows.Media;
#endif
namespace MapControl.Projections
{
/// <summary>
/// Spherical Orthographic Projection - AUTO2:42003.
/// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.148-150.
/// </summary>
public class Wgs84OrthographicProjection : ProjNetMapProjection
{
public Wgs84OrthographicProjection()
@ -25,5 +33,15 @@ namespace MapControl.Projections
CoordinateSystemWkt = string.Format(
CultureInfo.InvariantCulture, wktFormat, Center.Latitude, Center.Longitude);
}
public override Matrix RelativeScale(double latitude, double longitude)
{
var p = new AzimuthalProjection.ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude);
var h = p.CosC; // p.149 (20-5)
var scale = new Matrix(h, 0d, 0d, 1d, 0d, 0d);
scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI);
return scale;
}
}
}

View file

@ -1,7 +1,14 @@
using System.Globalization;
#if WPF
using System.Windows.Media;
#endif
namespace MapControl.Projections
{
/// <summary>
/// Spherical Stereographic Projection - AUTO2:97002.
/// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.157-160.
/// </summary>
public class Wgs84StereographicProjection : ProjNetMapProjection
{
public Wgs84StereographicProjection()
@ -25,5 +32,13 @@ namespace MapControl.Projections
CoordinateSystemWkt = string.Format(
CultureInfo.InvariantCulture, wktFormat, Center.Latitude, Center.Longitude);
}
public override Matrix RelativeScale(double latitude, double longitude)
{
var p = new AzimuthalProjection.ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude);
var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1
return new Matrix(k, 0d, 0d, k, 0d, 0d);
}
}
}