在 .NET 中使用 Proj4 进行坐标转换

• 下载电子表格结果 - 12.6 KB

引言

在我关于“转换地理坐标”（链接）的提示/技巧之后，我想讨论更高级的坐标转换主题。我将此提示/技巧标记为“高级”，因为它是面向高级用户受众的（代码相当简单）。如果您熟悉坐标转换和基准转换，那么此提示适合您。如果不是，您可以阅读：坐标系统、投影和基准转换，以及更深入的坐标转换。

背景

我是一名 GIS 用户很长时间了，一直对坐标转换及其背后的代码很感兴趣。很久以前，我甚至启动了我自己的“VBA 中的投影库”，但它已默默夭折。幸运的是，在这方面有一些非常好且免费的程序和 DLL 可用。我非常喜欢免费的 GIS 包 QGIS（仅代表我个人观点）。QGIS 和许多其他应用程序都使用名为 proj.4 的 C++ 库，该库已存在一段时间，经过了充分的测试，请参考 OSGeo Proj.4。据我所知，有两个 .NET 端口可用：Proj4Net 和 Dotspatial.Projections DLL，它是 DotSpatial 项目的一部分。我最初使用的是 Proj4Net，因为它是一个小巧、专用的 DLL，但发现它只支持几种标准投影后就放弃了。例如，如果您想将坐标从荷兰国家网格（称为 RDnew）转换为其他内容，您需要处理斜方位投影，而 Proj4Net 并不支持。然后，我转向了 DotSpatial.Projections 并成功实现了所有功能。该 DLL 有 19MB 大，但值得，它已针对 .NET4.0 编译。该 DLL 中有大量的方​​法，并且网上也有文档，因此我将不涵盖其所有功能。我只会展示它在 .NET 中的用法，并与荷兰皇家海军海道测量局开发的一个名为 PCTrans 的坐标转换包进行基准测试。该包已被证明对我来说非常准确；它是免费的，您可以 在此下载（但它在荷兰近海和近岸地区非常有用，您需要留下您的电子邮件地址）。

最后提醒一句：任何坐标转换或变换都是近似值。不存在精确的转换。一切都归结为转换误差是否足够好（小）以满足您的目标。如果您想精确，请使用您能获得的最好的 GPS 测量您的坐标，或者购买最好的 GPS，并且不要转换它们。

测试...

我将对 10 个点进行测试，源坐标为 ED50 基准的十进制度的纬度和经度，并将它们转换为荷兰陆地上使用的 RDnew 投影的东西向和南北向（X 和 Y）。ED50 基准使用 1924 年的国际椭球体（相当于 1909 年的“Hayford-Ellipsoid”）。RDnew 投影是基于 1841 年的 Bessel 椭球体的斜方位投影。上述转换意味着从 Int1924 椭球体到 Bessel1841 椭球体的基准转换。RDnew 投影仅在荷兰陆地部分有效（那里还有一个误差校正网格），但我将将其应用于近海区域，仅用于测试。下图显示了测试数据。

要引用 DotSpatial.Projections 库，我使用的是 SharpDevelop 的程序包管理器。您也可以通过 NuGet 库获取它。下图显示了程序包管理器。

编写程序是一件很容易的事（C# 代码在下面）

'install as package from https://nuget.net.cn/packages/DotSpatial.Projections/

Module Program
    
Sub Main()
    Console.WriteLine("VB.NET: Test coordinate conversion from ED50 to RD New")
    
    Dim x As Double() = {5.85000007,4.61923914,3.38732730,2.54133020,2.74920293, _
        3.38005007,5.10171662,6.03675746,6.83078962,6.77138124}
    
    Dim y As Double() = {50.84999997,51.51226471,51.43806765,51.89076392,54.36940251, _
        55.77197661,54.81463907,54.19125870,53.29330013,51.95008530}
    
    Dim z(x.Length-1) As Double
    
    For i As Integer = 0 To z.Length-1
        Console.WriteLine("input geographic ED50 p{0} = {1} {2}",i+1,x(i),y(i))
    Next
    
    'rewrite xy array for input into Proj4
    Dim xy(2*x.Length) As Double
    Dim ixy As Integer = 0
    For i As Integer = 0 To x.Length-1
        xy(ixy) = x(i)
        xy(ixy+1) = y(i)
        z(i) = 0
        ixy += 2
    Next
    
    Dim proj4_ed50 As String = "+proj=longlat +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +no_defs"
    Dim proj4_rdnew As String = "+proj=sterea " & _
       "+lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 " & _
       "+ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 " & _
       "+units=m +no_defs"
    
    Dim src As DotSpatial.Projections.ProjectionInfo = _
               DotSpatial.Projections.ProjectionInfo.FromProj4String(proj4_ed50)
    Dim trg As DotSpatial.Projections.ProjectionInfo = _
               DotSpatial.Projections.ProjectionInfo.FromProj4String(proj4_rdnew)
    
    DotSpatial.Projections.Reproject.ReprojectPoints(xy,z,src,trg,0,x.Length)

    ixy = 0
    For i As Integer = 0 To x.Length-1
        Console.WriteLine("output RD New {0} = {1} {2}",i+1,xy(ixy),xy(ixy+1))
        ixy += 2
    Next
    
    Console.Write("Press any key to continue . . . ")
    Console.ReadKey(True)
End Sub

End Module

首先，将 LatLon 加载到单独的数组中。然后，我将展示如何将各个数组重写为可以传递给 Proj4 的数组，Proj4 使用一种交错的 x 和 y 数组 {x1,y1,x2,y2....} 以及一个单独的 z 值数组（我将其设置为 0）。

设置投影很容易，只需声明一个 ProjectionInfo 类，并通过 string 定义源投影和目标投影。Proj4 允许您通过 4 种不同的方式定义投影：FromAuthorityCode、FromEpsgCode、FromEsriString 和 FromProj4String。有关这些方法的描述，请在此处查看：此处。我更喜欢 Proj4 string，因为它简洁、可读，而且更重要的是，它允许您调整投影参数，甚至可以轻松定义自己的自定义投影。这里的主要重点是 string 的 +towgs84 部分，它定义了基准转换。此标志后面的参数定义了当前基准如何转换为参考 WGS84 椭球体（“转换为 WGS84”）。定义这些后，就可以将坐标转换为不同的投影。Proj4 处理 3 参数转换（仅平移）和 7 参数（平移+旋转+缩放）转换。这 3 或 7 个数字的确切定义并非易事，具体取决于方法（Helmert、Bursa-Wolf、Molodensky 等）和旋转方向。网上有很多关于哪些参数最适合哪些区域的背景信息；只需搜索“proj4 3 or 7 parameters transformation”。很容易犯一个小小的拼写错误或使用错误的参数，导致坐标偏离您原以为的位置数百万英里……最佳实践是使用一个或一组您确信转换后的坐标是正确的测试点，并将其与您的结果进行比较。在我的例子中，我从 QGIS 中获取了 Proj4（包括 towgs84 参数）string，当您更改坐标参考系统时，QGIS 会记录这些信息。我的用于检查有效性的测试集就是这里展示的。从代码中可以看出，我使用了从 ED50 基准到 WGS84 的“简单”3 参数转换，以及从 Bessel1841 到 WGS84 的 7 参数转换。

DotSpatial.Projections.Reproject.ReprojectPoints(xy,z,src,trg,0,x.Length) 执行实际转换，其中 xy 和 z 是坐标数组，src 和 trg 是源投影和目标投影，0 和 x.Length 是数组的起始索引和长度。转换后的坐标将写入输入数组（或 VB.NET 中的 ByRef 参数）。

C# 代码在这里

using System;

namespace CPdemoInCSharp
{
    class Program
    {
        public static void Main(string[] args)
        {
        Console.WriteLine("C#: Test coordinate conversion from ED50 to RD New");

        double[] x = {5.85000007,4.61923914,3.3873273,2.5413302,
        2.74920293,3.38005007,5.10171662,6.03675746,6.83078962,6.77138124};
        double[] y = {50.84999997,51.51226471,51.43806765,
        51.89076392,54.36940251,55.77197661,54.81463907,54.1912587,53.29330013,51.9500853};
        double[] z = new double[x.Length];

        for (int i = 0; i <= z.Length - 1; i++) {
            Console.WriteLine("input geographic ED50 p{0} = {1} {2}", i+1, x[i], y[i]);
        }

        //rewrite xy array for input into Proj4
        double[] xy = new double[2 * x.Length];
        int ixy = 0;
        for (int i = 0; i <= x.Length - 1; i++) {
            xy[ixy] = x[i];
            xy[ixy + 1] = y[i];
            z[i] = 0;
            ixy += 2;
        }

        string proj4_ed50 = "+proj=longlat 
        +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +no_defs";
        string proj4_rdnew = "+proj=sterea +lat_0=52.15616055555555 
        +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel 
        +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs";

        DotSpatial.Projections.ProjectionInfo src = 
        	DotSpatial.Projections.ProjectionInfo.FromProj4String(proj4_ed50);
        DotSpatial.Projections.ProjectionInfo trg = 
        	DotSpatial.Projections.ProjectionInfo.FromProj4String(proj4_rdnew);

        DotSpatial.Projections.Reproject.ReprojectPoints(xy, z, src, trg, 0, x.Length);

        ixy = 0;
        for (int i = 0; i <= x.Length - 1; i++) {
            Console.WriteLine("output RD New p{0} = {1} {2}", i+1, xy[ixy], xy[ixy + 1]);
            ixy += 2;
        }

        Console.Write("Press any key to continue . . . ");
        Console.ReadKey(true);
        }
    }
}

结果...

转换的结果和比较列在下表（也附在本技巧的电子表格中）。

可以看出，与 PCTrans 的输出相比，结果非常接近。当 PCTrans 将坐标转换为 RDnew 时，它会警告输入坐标是否超出有效性和校正网格可用区域。在输入坐标位于有效区域内时，绝对差异最大为 18 厘米。请注意，在此有效区域内，PCTrans 应用了额外的校正。在有效区域外（此时 PCTrans不应用任何校正），差异小于 1 厘米。对于我的工作来说，这种精度已经足够了。如果您需要毫米级的精度，这些结果可能不会让您满意；但请记住，转换本身就是一种近似，并且永远不会 100% 准确。

希望这个技巧对您有用...

历史

• 2016 年 1 月 18 日：初稿