diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs
index fc08d8b2..193795b4 100644
--- a/MapControl/Shared/AzimuthalEquidistantProjection.cs
+++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs
@@ -26,7 +26,7 @@ namespace MapControl
CrsId = crsId;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var p = GetProjectedPoint(latitude, longitude);
diff --git a/MapControl/Shared/EquirectangularProjection.cs b/MapControl/Shared/EquirectangularProjection.cs
index b34270b2..6e4fea1b 100644
--- a/MapControl/Shared/EquirectangularProjection.cs
+++ b/MapControl/Shared/EquirectangularProjection.cs
@@ -28,7 +28,7 @@ namespace MapControl
CrsId = crsId;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
return new Matrix(1d / Math.Cos(latitude * Math.PI / 180d), 0d, 0d, 1d, 0d, 0d);
}
diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs
index f03c15c0..59d57c7c 100644
--- a/MapControl/Shared/GnomonicProjection.cs
+++ b/MapControl/Shared/GnomonicProjection.cs
@@ -26,7 +26,7 @@ namespace MapControl
CrsId = crsId;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var p = GetProjectedPoint(latitude, longitude);
var k = 1d / p.CosC; // p.165 (22-3)
diff --git a/MapControl/Shared/MapBase.cs b/MapControl/Shared/MapBase.cs
index c45be0ca..321936d9 100644
--- a/MapControl/Shared/MapBase.cs
+++ b/MapControl/Shared/MapBase.cs
@@ -183,7 +183,7 @@ namespace MapControl
///
public Matrix GetMapToViewTransform(double latitude, double longitude)
{
- var transform = MapProjection.RelativeScale(latitude, longitude);
+ var transform = MapProjection.RelativeTransform(latitude, longitude);
transform.Scale(ViewTransform.Scale, ViewTransform.Scale);
transform.Rotate(ViewTransform.Rotation);
diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs
index fcdb89f2..3aaa3fa8 100644
--- a/MapControl/Shared/MapProjection.cs
+++ b/MapControl/Shared/MapProjection.cs
@@ -103,10 +103,10 @@ namespace MapControl
protected virtual void CenterChanged() { }
///
- /// Gets the relative scale at the specified geographic coordinates.
+ /// Gets the relative transform at the specified geographic coordinates.
/// The returned Matrix represents the local distortion of the map projection.
///
- public abstract Matrix RelativeScale(double latitude, double longitude);
+ public abstract Matrix RelativeTransform(double latitude, double longitude);
///
/// Transforms geographic coordinates to a Point in projected map coordinates.
@@ -121,9 +121,9 @@ namespace MapControl
public abstract Location MapToLocation(double x, double y);
///
- /// Gets the relative map scale at the specified geographic Location.
+ /// Gets the relative transform at the specified geographic Location.
///
- public Matrix RelativeScale(Location location) => RelativeScale(location.Latitude, location.Longitude);
+ public Matrix RelativeTransform(Location location) => RelativeTransform(location.Latitude, location.Longitude);
///
/// Transforms a Location in geographic coordinates to a Point in projected map coordinates.
diff --git a/MapControl/Shared/Matrix.cs b/MapControl/Shared/Matrix.cs
index c3ffeacf..b54b7bc9 100644
--- a/MapControl/Shared/Matrix.cs
+++ b/MapControl/Shared/Matrix.cs
@@ -52,8 +52,6 @@ namespace MapControl
public void Scale(double scaleX, double scaleY)
{
- // equivalent to Multiply(new Matrix(scaleX, 0, 0, scaleY, 0d, 0d));
- //
SetMatrix(
M11 * scaleX,
M12 * scaleY,
@@ -72,8 +70,6 @@ namespace MapControl
var cos = Math.Cos(angle);
var sin = Math.Sin(angle);
- // equivalent to Multiply(new Matrix(cos, sin, -sin, cos, 0d, 0d));
- //
SetMatrix(
M11 * cos - M12 * sin,
M11 * sin + M12 * cos,
@@ -99,15 +95,15 @@ namespace MapControl
invDet * (M12 * OffsetX - M11 * OffsetY));
}
- public void Multiply(Matrix m)
+ public static Matrix Multiply(Matrix m1, Matrix m2)
{
- SetMatrix(
- M11 * m.M11 + M12 * m.M21,
- M11 * m.M12 + M12 * m.M22,
- M21 * m.M11 + M22 * m.M21,
- M21 * m.M12 + M22 * m.M22,
- OffsetX * m.M11 + OffsetY * m.M21 + m.OffsetX,
- OffsetX * m.M12 + OffsetY * m.M22 + m.OffsetY);
+ return new Matrix(
+ m1.M11 * m2.M11 + m1.M12 * m2.M21,
+ m1.M11 * m2.M12 + m1.M12 * m2.M22,
+ m1.M21 * m2.M11 + m1.M22 * m2.M21,
+ m1.M21 * m2.M12 + m1.M22 * m2.M22,
+ m1.OffsetX * m2.M11 + m1.OffsetY * m2.M21 + m2.OffsetX,
+ m1.OffsetX * m2.M12 + m1.OffsetY * m2.M22 + m2.OffsetY);
}
private void SetMatrix(double m11, double m12, double m21, double m22, double offsetX, double offsetY)
diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs
index a09a7059..263576a0 100644
--- a/MapControl/Shared/OrthographicProjection.cs
+++ b/MapControl/Shared/OrthographicProjection.cs
@@ -26,7 +26,7 @@ namespace MapControl
CrsId = crsId;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var p = GetProjectedPoint(latitude, longitude);
diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs
index 71276c1d..787b5d73 100644
--- a/MapControl/Shared/PolarStereographicProjection.cs
+++ b/MapControl/Shared/PolarStereographicProjection.cs
@@ -39,7 +39,7 @@ namespace MapControl
return r / m; // p.161 (21-32)
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var k = RelativeScale(Hemisphere, Flattening, latitude);
diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs
index ce7cab88..9b891793 100644
--- a/MapControl/Shared/StereographicProjection.cs
+++ b/MapControl/Shared/StereographicProjection.cs
@@ -26,7 +26,7 @@ namespace MapControl
CrsId = crsId;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var p = GetProjectedPoint(latitude, longitude);
var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1
diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs
index 5ab6a681..3ab4da9c 100644
--- a/MapControl/Shared/TransverseMercatorProjection.cs
+++ b/MapControl/Shared/TransverseMercatorProjection.cs
@@ -60,24 +60,21 @@ namespace MapControl
Flattening = Wgs84Flattening;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
- // h = k/k0 for relative scale
- var h = 1d;
-
-#if TM_RELATIVE_SCALE // relative scale is usually < 1.001 and hence neglectable
// φ
var phi = latitude * Math.PI / 180d;
var sinPhi = Math.Sin(phi);
- // t
- var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi));
// λ - λ0
var dLambda = (longitude - CentralMeridian) * Math.PI / 180d;
var cosLambda = Math.Cos(dLambda);
+ // t
+ var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi));
+ var u = Math.Sqrt(1d + t * t);
// ξ'
var xi_ = Math.Atan2(t, cosLambda);
// η'
- var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t));
+ var eta_ = Atanh(Math.Sin(dLambda) / u);
// σ
var sigma = 1 +
2d * a1 * Math.Cos(2d * xi_) * Math.Cosh(2d * eta_) +
@@ -90,10 +87,18 @@ namespace MapControl
6d * a3 * Math.Sin(6d * xi_) * Math.Sinh(6d * eta_);
var n = Flattening / (2d - Flattening);
- var u = (1d - n) / (1d + n) * Math.Tan(phi);
- h = f1 * Math.Sqrt((1d + u * u) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda));
-#endif
- return new Matrix(h, 0d, 0d, h, 0d, 0d);
+ var m = (1d - n) / (1d + n) * Math.Tan(phi);
+ // h = k/k0 for relative scale
+ var h = f1 * Math.Sqrt((1d + m * m) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda));
+
+ var tanLambda = Math.Tan(dLambda);
+ // γ, meridian convergence angle
+ var gamma = Math.Atan2(tau * u + sigma * t * tanLambda, sigma * u - tau * t * tanLambda);
+
+ var transform = new Matrix(h, 0d, 0d, h, 0d, 0d);
+ transform.Rotate(-gamma * 180d / Math.PI);
+
+ return transform;
}
public override Point? LocationToMap(double latitude, double longitude)
diff --git a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
index c015dac0..e9ab561d 100644
--- a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
+++ b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
@@ -45,23 +45,26 @@ namespace MapControl
e6 * 35d / 3072d * Math.Sin(6d * phi)); // (3-21)
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
// h = k/k0 for relative scale
var h = 1d;
+ // γ, meridian convergence angle
+ var gamma = 0d;
-#if TM_RELATIVE_SCALE // relative scale is usually < 1.001 and hence neglectable
if (latitude > -90d && latitude < 90d)
{
var phi = latitude * Math.PI / 180d;
+ var sinPhi = Math.Sin(phi);
var cosPhi = Math.Cos(phi);
- var tanPhi = Math.Tan(phi);
+ var tanPhi = sinPhi / cosPhi;
+ var dLambda = (longitude - CentralMeridian) * Math.PI / 180d;
var e2 = (2d - Flattening) * Flattening;
var e_2 = e2 / (1d - e2); // (8-12)
var T = tanPhi * tanPhi; // (8-13)
var C = e_2 * cosPhi * cosPhi; // (8-14)
- var A = (longitude - CentralMeridian) * Math.PI / 180d * cosPhi; // (8-15)
+ var A = dLambda * cosPhi; // (8-15)
var A2 = A * A;
var A4 = A2 * A2;
var A6 = A2 * A4;
@@ -69,9 +72,14 @@ namespace MapControl
h = 1d + (1d + C) * A2 / 2d +
(5d - 4d * T + 42d * C + 13d * C * C - 28d * e_2) * A4 / 24d +
(61d - 148d * T + 16 * T * T) * A6 / 720d; // (8-11)
+
+ gamma = Math.Atan(Math.Tan(dLambda) * sinPhi);
}
-#endif
- return new Matrix(h, 0d, 0d, h, 0d, 0d);
+
+ var transform = new Matrix(h, 0d, 0d, h, 0d, 0d);
+ transform.Rotate(-gamma * 180d / Math.PI);
+
+ return transform;
}
public override Point? LocationToMap(double latitude, double longitude)
diff --git a/MapControl/Shared/WebMercatorProjection.cs b/MapControl/Shared/WebMercatorProjection.cs
index 528a0200..5c5149f7 100644
--- a/MapControl/Shared/WebMercatorProjection.cs
+++ b/MapControl/Shared/WebMercatorProjection.cs
@@ -27,7 +27,7 @@ namespace MapControl
CrsId = crsId;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var k = 1d / Math.Cos(latitude * Math.PI / 180d); // p.44 (7-3)
diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs
index 7f06dee9..8a131872 100644
--- a/MapControl/Shared/WorldMercatorProjection.cs
+++ b/MapControl/Shared/WorldMercatorProjection.cs
@@ -29,7 +29,7 @@ namespace MapControl
CrsId = crsId;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var k = RelativeScale(latitude);
diff --git a/MapProjections/Shared/ProjNetMapProjection.cs b/MapProjections/Shared/ProjNetMapProjection.cs
index f18c520c..02c57c7b 100644
--- a/MapProjections/Shared/ProjNetMapProjection.cs
+++ b/MapProjections/Shared/ProjNetMapProjection.cs
@@ -71,7 +71,7 @@ namespace MapControl.Projections
public MathTransform MapToLocationTransform { get; private set; }
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
return new Matrix(1d, 0d, 0d, 1d, 0d, 0d);
}
diff --git a/MapProjections/Shared/WebMercatorProjection.cs b/MapProjections/Shared/WebMercatorProjection.cs
index 36aa9a68..9c32b71d 100644
--- a/MapProjections/Shared/WebMercatorProjection.cs
+++ b/MapProjections/Shared/WebMercatorProjection.cs
@@ -17,7 +17,7 @@ namespace MapControl.Projections
CoordinateSystem = ProjectedCoordinateSystem.WebMercator;
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var k = 1d / Math.Cos(latitude * Math.PI / 180d); // p.44 (7-3)
diff --git a/MapProjections/Shared/Wgs84OrthographicProjection.cs b/MapProjections/Shared/Wgs84OrthographicProjection.cs
index fe89dcd9..400f6f60 100644
--- a/MapProjections/Shared/Wgs84OrthographicProjection.cs
+++ b/MapProjections/Shared/Wgs84OrthographicProjection.cs
@@ -34,7 +34,7 @@ namespace MapControl.Projections
CultureInfo.InvariantCulture, wktFormat, Center.Latitude, Center.Longitude);
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var p = new AzimuthalProjection.ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude);
diff --git a/MapProjections/Shared/Wgs84StereographicProjection.cs b/MapProjections/Shared/Wgs84StereographicProjection.cs
index 3da4d8df..c697007f 100644
--- a/MapProjections/Shared/Wgs84StereographicProjection.cs
+++ b/MapProjections/Shared/Wgs84StereographicProjection.cs
@@ -34,7 +34,7 @@ namespace MapControl.Projections
CultureInfo.InvariantCulture, wktFormat, Center.Latitude, Center.Longitude);
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(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
diff --git a/MapProjections/Shared/Wgs84UpsProjections.cs b/MapProjections/Shared/Wgs84UpsProjections.cs
index 36d14f0b..fdb6acb5 100644
--- a/MapProjections/Shared/Wgs84UpsProjections.cs
+++ b/MapProjections/Shared/Wgs84UpsProjections.cs
@@ -21,7 +21,7 @@ namespace MapControl.Projections
"AUTHORITY[\"EPSG\",\"32661\"]]";
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var k = PolarStereographicProjection.RelativeScale(Hemisphere.North, Wgs84Flattening, latitude);
@@ -46,7 +46,7 @@ namespace MapControl.Projections
"AUTHORITY[\"EPSG\",\"32761\"]]";
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var k = PolarStereographicProjection.RelativeScale(Hemisphere.South, Wgs84Flattening, latitude);
diff --git a/MapProjections/Shared/WorldMercatorProjection.cs b/MapProjections/Shared/WorldMercatorProjection.cs
index 8a6b6352..63c9708c 100644
--- a/MapProjections/Shared/WorldMercatorProjection.cs
+++ b/MapProjections/Shared/WorldMercatorProjection.cs
@@ -27,7 +27,7 @@ namespace MapControl.Projections
"AUTHORITY[\"EPSG\",\"3395\"]]";
}
- public override Matrix RelativeScale(double latitude, double longitude)
+ public override Matrix RelativeTransform(double latitude, double longitude)
{
var k = MapControl.WorldMercatorProjection.RelativeScale(latitude);