WGS84 转瑞典国家网格 (RT90) 投影类
一篇关于如何将平面坐标从 WGS 84 椭球体转换到瑞典国家网格 (Riket koordinatsystem 1990) 的文章。
引言
这篇简短的文章描述了如何将 NMEA 字符串（即来自 GPS 设备）转换为瑞典国家网格 (RT90)。 我不会详细介绍 SerialPort 通信或 NMEA 解释问题，因为 Jon Person 在以下文章中对此进行了很好的描述：
背景
当我在工作时编写我的物理应用程序（包括与 GPS 设备通信）时，我遇到了上面列出的文章。 Jon Person 的 NMEA-Interpreter 文章特别有帮助（谢谢 Jon），但自然地，它们都没有涵盖如何将 NMEA 字符串（基于 WGS 84 椭球体）转换为瑞典国家网格（Rikets Koordinatsystem 1990 或简称 RT90），这似乎是瑞典大多数地图使用的。 我还在Lantmäteriet找到了很多有用的信息，其中提供了投影/变换公式。 所需要的只是实现转换，最好是用 C#。 所以，这就是。
使用代码
该代码只是一个简单的类，它读取 NMEA 字符串，例如 Jon Person 的 NMEA-Interpreter 给出的那些字符串，并将它们转换为 RT90 笛卡尔坐标 (x, y)。 我还包括了对瑞典 6 个不同地区的支持，范围从东部的吕勒奥到西部的哥德堡，正如 Lantmäteriet 建议的那样。 这些区域被称为 5 gon O、2.5 gon O、0 gon V、2.5 gon V、5 gon V 和 7.5 gon V。主要（在瑞典意义上）的城市从东到西依次是：吕勒奥、于默奥、斯德哥尔摩、厄勒布鲁、马尔默和哥德堡。（注意：'V' 和 'O' 代替 'W' 和 'E'，因为瑞典语中西/东的词以 'V'/'O' 开头......）
该类被称为 GaussKrugerProjection，因为它使用的是保角投影方法，即“保留感兴趣区域的原始形状但不保留面积或距离的投影”。 它的使用方法很简单，如下所示：
// Create a new Transform instance, using 0 gon V (Stockholm):
GaussKrugerProjection myTransformInStockholm = new GaussKrugerProjection("0V");
或者，如果您在哥德堡 (Gothenburg) 附近：
// Create a new Transform instance, using 7.5 gon V:
GaussKrugerProjection myTransformInGothenburg = new
GaussKrugerProjection("7.5V");
如果我们收到了来自例如上面提到的 NMEA-Interpreter 的 NMEA 字符串，我们只需使用 GaussKrugerProjection 类中的 GetRT90 方法将其发送到投影：
public int[] GetRT90(string lat, string lon)
...使用上面创建的投影实例，并假设我们有一个 NMEAInterpreter
//Get string from NMEA Interpreter
private void NMEAInterpreter_PositionReceived(string lat, string lon)
{
// Get x and y coordinates in RT90 0 gon V
int[] XandY = myTransformInStockholm.GetRT90(lat, lon);
//Where x is at index 0 and y at index 1 in the int[]
int x = XandY[0];
int y = XandY[1];
// OK, here we have the RT90 0 gon V coordinates
Console.WriteLine("x = {0}, y = {1}", x.ToString(), y.ToString());
}
该类是为以下格式的 NMEA 字符串编写的：
// Latitude example
string lat = "20°18.2274\"E";
其中 20 是小时，18 分钟，22.74 秒，E 代表半球。 确保您的 NMEA 字符串以相同的方式格式化。
GaussKrugerProjection 类基本上分为两个部分：
- NMEA 字符串到弧度解析器（类似于上面文章 (3) 中的解析器）
- 投影计算器。
解析器看起来像这样：
// isLong is true for Longitudes, otherwise false
private double GetLatLong(string LatorLong, bool isLong)
{
//Get Hours (up to the '°')
double deciLatLon = Convert.ToDouble(
LatorLong.Substring(0, LatorLong.IndexOf("°"))
);
//Remove it once we've used it
LatorLong = LatorLong.Substring(LatorLong.IndexOf("°") + 1);
//Get Minutes (up to the '.') and divide by Minutes/Hour
deciLatLon +=
(Convert.ToDouble(LatorLong.Substring(
0, LatorLong.IndexOf(".")), enUSCultureInfo)
) / 60.0;
//Remove it once we've used it
LatorLong = LatorLong.Substring(LatorLong.IndexOf(".") + 1);
//Get Seconds (up to the '"') and divide by Seconds/Hour
string sec = LatorLong.Substring(0, LatorLong.IndexOf("\""));
// Insert a dot after the first two digits in seconds
deciLatLon +=
(
Convert.ToDouble(sec.Insert(2, "."), enUSCultureInfo)
) / 3600.0;
//Get the Hemisphere
LatorLong = LatorLong.Substring(
LatorLong.IndexOf("\"") + 1);
if (isLong && LatorLong == "S" ||
!isLong && LatorLong == "W")
{
// Set us right if we're not in the
// West || Northern hemisphere
deciLatLon = -deciLatLon;
}
//And return (as radians)
return deciLatLon * (Math.PI / 180.0);
}
投影计算器可以分为两个部分：
- 声明变量/预投影计算和
- 执行投影。
首先：这是变量声明
#region Fields
//GRS 80 Ellipsoid Characteristics:
//Semi Major axis
private static double a = 6378137.0;
//Flattening
private static double f = (1.0 / 298.2572221010);
//RT90 0 gon V 0:-15 fields (Use around Stockholm)
// Centrum meridian
private static string CM_0V = "18°03.2268\"E";
// Scale factor
private static double k0_0V = 1.000005400000;
// False North
private static double FN_0V = -668.844;
// False East
private static double FE_0V = 1500083.521;
//Variables
private string CM;
private double k0;
private double FN;
private double FE;
private double lat; // Geodetic latitude
private double lon; // Geodetic longitude
//Gauss-Krüger Projection variables
private double A, B, C, D, Beta1, Beta2, Beta3, Beta4,
e2, n, aHat;
//RT90-coordinates
private double x;
private double y;
//Make it international...(for numers in NMEA-sentences)
private static CultureInfo enUSCultureInfo =
new CultureInfo("en-US");
#endregion
对于所有 6 个不同的 gon，都存在类似的 Fields。 然后构造函数使用它们：
/// <summary>
/// Good for use in the Stockholm-region.
/// </summary>
public GaussKrugerProjection()
{
// USE 2.5 V 0:-15
CM = CM_25V;
k0 = k0_25V;
FN = FN_25V;
FE = FE_25V;
this.Initialize();
}
private void Initialize()
{
e2 = f * (2.0 - f);
n = f / (2.0 - f);
aHat = (a / (1.0 + n)) * (1.0 +
(0.25 * Math.Pow(n, 2)) +
((1.0 / 64.0) * Math.Pow(n, 4)));
A = e2;
B = (1.0 / 6.0) * (5.0 * Math.Pow(A, 2) -
6.0 * Math.Pow(A, 3));
C = (1.0 / 120.0) * (104.0 * Math.Pow(A, 3) -
45.0 * Math.Pow(A, 4));
D = (1.0 / 1260.0) * (1237.0 * Math.Pow(A, 4));
Beta1 = (0.5 * n) - ((2.0 / 3.0) * Math.Pow(n, 2)) +
((5.0 / 16.0) * Math.Pow(n, 3)) +
((41.0 / 180.0) * Math.Pow(n, 4));
Beta2 = ((13.0 / 48.0) * Math.Pow(n, 2)) -
((3.0 / 5.0) * Math.Pow(n, 3)) +
((557.0 / 1440.0) * Math.Pow(n, 4));
Beta3 = ((61.0 / 240.0) * Math.Pow(n, 3)) -
((103.0 / 140.0) * Math.Pow(n, 4));
Beta4 = ((49561.0 / 161280.0) * Math.Pow(n, 4));
}
最后是实际投影
private void CalcGaussKrugerProjection(double lat, double lon)
{
// Compute the Conformal Latitude
double phiStar = lat - (Math.Sin(lat) * Math.Cos(lat) * (
A +
B * Math.Pow(Math.Sin(lat), 2) +
C * Math.Pow(Math.Sin(lat), 4) +
D * Math.Pow(Math.Sin(lat), 6)));
// Difference in longitude
double dLon = lon - GetLatLong(CM, true);
// Get Angles:
double chi = Math.Atan(Math.Tan(phiStar) / Math.Cos(dLon));
//Since Atanh isn't represented in the Math-class
//we'll use a simplification that holds for real z < 1
//Ref:
//https://mathworld.net.cn/InverseHyperbolicTangent.html
double z = Math.Cos(phiStar) * Math.Sin(dLon);
double eta = 0.5 * Math.Log((1.0 + z) / (1.0 - z));
// OK , we're finally ready to calculate the
// cartesian coordinates in RT90
x = k0 * aHat * (chi +
Beta1 * Math.Sin(2.0 * chi) * Math.Cosh(2.0 * eta) +
Beta2 * Math.Sin(4.0 * chi) * Math.Cosh(4.0 * eta) +
Beta3 * Math.Sin(6.0 * chi) * Math.Cosh(6.0 * eta) +
Beta4 * Math.Sin(8.0 * chi) * Math.Cosh(8.0 * eta)) +
FN;
y = k0 * aHat * (eta +
Beta1 * Math.Cos(2.0 * chi) * Math.Sinh(2.0 * eta) +
Beta2 * Math.Cos(4.0 * chi) * Math.Sinh(4.0 * eta) +
Beta3 * Math.Cos(6.0 * chi) * Math.Sinh(6.0 * eta) +
Beta4 * Math.Cos(8.0 * chi) * Math.Sinh(8.0 * eta)) +
FE;
}
如果您对数学细节感兴趣，请参阅：此页面或此 PDF了解更多信息（第一个是瑞典语）。
历史
这是第一个版本。
- 欢呼！