大地变换

PROJ可以做任何事情,从最简单的投影到跨越许多参考帧的非常复杂的变换。虽然PROJ最初是作为地图投影工具开发的,但随着时间的推移,它已经演变成一个强大的通用坐标转换引擎,使其能够在大地高精度水平上进行大规模地图投影和坐标转换。本章深入研究了如何进行复杂程度不同的大地测量变换的细节。

在PROJ中,存在两种大地坐标变换框架,即 PROJ 4.x/5.x / cs2cs / pj_transform() 框架和 变压管道 框架。第一个是原始的,有限的,用于在PROJ中进行大地测量转换的框架,后者是一个更新的补充,旨在成为一个更完整的转换框架。以下各节对这两种情况进行了说明。正文的很大一部分是基于 [EversKnudsen2017] .

在描述这两个框架的细节之前,让我们首先指出,大地坐标变换的大多数情况可以表示为一系列基本运算,一个运算的输出就是下一个运算的输入。E、 g.从UTM 32区基准ED50到UTM 32区基准ETRS89时,在最简单的情况下,必须经过5个步骤:

  1. 将UTM坐标反向投影到地理坐标

  2. 将地理坐标转换为三维笛卡尔地心坐标

  3. 应用从ED50到ETRS89的Helmert变换

  4. 从笛卡尔坐标转换回地理坐标

  5. 最后将地理坐标投影到UTM zone 32平面坐标。

变压管道

上述步骤与unixshell样式管道之间的同源性是显而易见的。这就是转换管道框架背后的主要建筑灵感。流水线框架是通过使用一个特殊的“投影”来实现的,该投影将一系列基本操作作为用户提供的参数,并将这些操作串在一起以实现所需的全部转换。此外,许多基本的大地测量操作,包括Helmert变换、一般高阶多项式变换和Molodensky变换,都可以作为管道框架的一部分。在预期即将到来的支持全时变转换,我们还介绍了一个4D时空数据类型,和一个编程接口(API)来处理这个。

莫洛登斯基变换直接从一个基准中的大地坐标转换为另一个基准中的大地坐标,而赫尔默特变换(通常更精确)则从三维笛卡尔坐标转换为三维笛卡尔坐标。因此,当使用Helmert变换时,通常需要进行从大地坐标到笛卡尔坐标的初始转换,以及相反方向的最终转换,以获得所需的结果。幸运的是,这种三步复合变换有一个吸引人的特点,即每一步只依赖于前一步的输出。因此,我们可以通过将三个步骤(大地坐标到笛卡尔→Helmert→笛卡尔坐标到大地坐标)的输出和输入连接起来,以流水线方式构建大地坐标到大地坐标的Helmert变换。管道驱动程序使这种链式转换成为可能。实现非常紧凑,只包含一个伪投影,称为 pipeline ,它将基本投影的字符串作为其参数(注意:“投影”是用于任何类型转换的一个有点误导性的PROJ术语)。管道伪投影由许多基本变换补充,总的来说提供了一个框架,用于为广泛的大地测量任务构建高精度解决方案。

作为第一个例子,让我们来看看标志性的 大地测量→笛卡尔→赫尔默特→大地测量 案例(介绍中示例中的步骤2到4)。在项目中,它可以实现为

proj=pipeline
step proj=cart ellps=intl
step proj=helmert convention=coordinate_frame
     x=-81.0703  y=-89.3603  z=-115.7526
    rx=-0.48488 ry=-0.02436 rz=-0.41321  s=-0.540645
step proj=cart inv ellps=GRS80

管道可以在两端扩展,以适应输入和输出所需的任何坐标类型:在下面的示例中,我们将不推荐使用的丹麦系统45(在原始定义网络中具有一定张力的2D系统)转换为UTM区域33,ETRS89。使用多项式变换降低张力(初始值=./s45b。。。步骤s45b.pol是一个包含多项式系数的文件),将S45坐标转换为技术坐标系(TC32),定义为表示“UTM区域32坐标,因为它们将查看ED50和ETRS89之间的Helmert转换是否完美”。然后使用UTM逆投影将TC32坐标转换回大地坐标(ED50),再转换回笛卡尔坐标(ED50),然后使用相关的Helmert变换转换回笛卡尔坐标(ETRS89),再转换回大地坐标(ETRS89),最后投影到UTM 33区ETRS89系统。总而言之,这是一个6步流水线,实现了一个厘米级精度的转换,该转换来自一个不推荐使用的分米级张力系统。

proj=pipeline
step init=./s45b.pol:s45b_tc32
step proj=utm inv ellps=intl zone=32
step proj=cart ellps=intl
step proj=helmert convention=coordinate_frame
      x=-81.0703  y=-89.3603  z=-115.7526
     rx=-0.48488 ry=-0.02436 rz=-0.41321 s=-0.540645
step proj=cart inv ellps=GRS80
step proj=utm ellps=GRS80 zone=33

有了管道框架,时空转换成为可能。这可以通过利用PROJ中的时间维度来实现,该维度允许4D坐标(三个空间分量和一个时间分量)通过转换管道。下面是从f2000到f93的转换示例。时间分量在输入数据中表示为GPS周,但14参数Helmert变换期望时间单位为小数年。因此,管道中的第一步是unitconvert伪投影,它确保将正确的单位传递给Helmert变换。Helmert变换的大部分参数取自 [Altamimi2002] ,除了转变的时代。管道中的最后一步是将坐标时间戳转换回GPS周。

proj=pipeline
step proj=unitconvert t_in=gps_week t_out=decimalyear
step proj=helmert convention=coordinate_frame
     x=0.0127 y=0.0065 z=-0.0209 s=0.00195
     rx=0.00039 ry=-0.00080 rz=0.00114
     dx=-0.0029 dy=-0.0002 dz=-0.0006 ds=0.00001
     drx=0.00011 dry=0.00019 drz=-0.00007
     t_epoch=1988.0
step proj=unitconvert t_in=decimalyear t_out=gps_week

项目4.x/5.x范例

参数

描述

+基准面

基准名称(参见 proj -ld

+土工格栅

用于垂直基准转换的GTX网格文件的文件名

+nadgrids公司

用于基准转换的NTv2网格文件的文件名

+拖车84

3或7项基准转换参数

+to_meter

将地图单位转换为1.0m的乘数

+vto_meter

垂直转换为米

警告

本节记录了PROJ 6.x中PROJ 4.x和5.x的行为, cs2cs 已重新加工以使用 proj_create_crs_to_crs() 内部,带有 后期装订 功能,因此不再局限于使用WGS84作为轴心(也称为 早期装订 方法)。什么时候? cs2cs PROJ.4扩展字符串用于描述CRS,包括 +towgs84+nadgrids+geoidgrids ,它通常会给出与早期项目版本相同的结果。与一起使用时机构:代码CRS描述,它可能返回不同的结果。

这个 CS2CS proj4和proj5中的框架提供了 管道 框架。在这一过程中,两个坐标变换作为一个坐标轴进行。也就是说,输入坐标转换为WGS84大地坐标,然后通过使用Helmert变换、基准移动网格或两者的组合,从WGS84坐标转换为指定的输出坐标参考系。基准位移可以用带有参数的项目串来描述 +towgs84+nadgrids+geoidgrids . 三者都存在逆变换,如果在输入项目字符串中指定,则应用逆变换。最常见的是 +towgs84 ,用于定义从输入参考坐标系到WGS84的3或7参数Helmert偏移。WGS84的具体实现没有具体说明,因此在转换的这一步中引入了相当多的不确定性。使用+nadgrids参数,可以应用由校正网格中的插值导出的非线性平面校正。最初,这是一种在北美基准NAD27和NAD83之间转换坐标的方法,但校正可以应用于存在校正网格的任何基准。水平网格偏移的逆变换是“哑的”,即逐字应用校正网格而不考虑逆运算是非线性的。类似于水平网格校正, +geoidgrids 可用于在垂直分量中执行栅格校正。两种栅格校正方法都允许在同一变换中包含多个栅格

与之相对的是 改造管道 框架,使用 CS2cs 项目4和项目5中的框架被表示为两个单独的项目字符串。一条投影线 to WGS84和One from WGS84。它们一起形成从源坐标参考系到目的坐标参考系的映射。当与 cs2cs 源CRS和目标CRS由特殊的 +to 参数。

用下面的示例将RSWGS887转换为WGG4 +towgs84 参数。

cs2cs +proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62
    +to +proj=latlong +datum=WGS84
20 35
20d0'5.467"E    35d0'9.575"N 0.000

使用PROJ 6,您可以简单地使用以下内容:

备注

对于PROJ 6,EPSG地理坐标参考系的坐标顺序是纬度第一,经度第二。

cs2cs "GGRS87" "WGS 84"
35 20
35d0'9.575"N    20d0'5.467"E 0.000

cs2cs EPSG:4121 EPSG:4326
35 20
35d0'9.575"N    20d0'5.467"E 0.000

EPSG数据库提供了使用近似的7参数转换从WGS72转换到WGS84的示例。

cs2cs +proj=latlong +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.219 \
    +to +proj=latlong +datum=WGS84
4 55
4d0'0.554"E     55d0'0.09"N 0.000

对于proj6,您可以简单地使用以下命令(注意纬度和经度的相反顺序)

cs2cs "WGS 72" "WGS 84"
55 4
55d0'0.09"N 4d0'0.554"E 0.000

cs2cs EPSG:4322 EPSG:4326
55 4
55d0'0.09"N 4d0'0.554"E 0.000

基于栅格的基准调整

在许多地方(特别是北美和澳大利亚),国家大地测量组织提供网格偏移文件,用于在不同基准之间进行转换,如NAD27到NAD83。这些网格偏移文件包括要应用于每个网格位置的偏移。实际上,通常基于包含四个网格点的网格之间的插值来计算网格偏移。

Proj支持使用栅格文件在各种参考系之间移动。格网平移表格格式为CTable、NTv1(旧的加拿大格式)和NTv2 (.gsb -新的加拿大和澳大利亚格式)。

本节中的文本基于 CS2CS 框架。网格转移也可能与 管道 框架。两者的主要区别在于 CS2CS 框架仅限于到WGS84的网格映射,而 变压管道 只要存在网格,就可以在任意两个参考帧之间执行网格移动。

将栅格移动用于 cs2cs 是使用 +nadgrids 坐标系定义中的关键字。例如:

% cs2cs +proj=latlong +ellps=clrk66 +nadgrids=ntv1_can.dat \
    +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 50
EOF
111d0'2.952"W   50d0'0.111"N 0.000

在本例中, /usr/local/share/proj/ntv1_can.dat 加载栅格移位文件,并使用该文件获取选定点的栅格移动值。

可以列出多个网格移位文件,在这种情况下,将依次尝试每个文件,直到找到一个包含要转换的点。

cs2cs +proj=latlong +ellps=clrk66 \
          +nadgrids=conus,alaska,hawaii,stgeorge,stlrnc,stpaul \
    +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 44
EOF
111d0'2.788"W   43d59'59.725"N 0.000

跳过丢失的网格

特殊前缀 @ 可以作为栅格的前缀以使其成为可选的。如果没有找到,搜索将继续到下一个网格。正常情况下,任何找不到的网格都会导致错误。例如,以下代码将使用 ntv2_0.gsb 文件,否则它将回退到使用 ntv1_can.dat 文件。

cs2cs +proj=latlong +ellps=clrk66 +nadgrids=@ntv2_0.gsb,ntv1_can.dat \
    +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 50
EOF
111d0'3.006"W   50d0'0.103"N 0.000

空网格

一份特价套餐 null 网格平移文件与Proj一起分发。这个文件提供了整个世界的零位移值。如果您想要将零位移应用于所有其他网格的有效区域之外的点,则可将其列在nadgrids文件列表的末尾。通常,如果找不到包含要转换的点的网格,则会出现错误。

cs2cs +proj=latlong +ellps=clrk66 +nadgrids=conus,null \
    +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 45
EOF
111d0'3.006"W   50d0'0.103"N 0.000

cs2cs +proj=latlong +ellps=clrk66 +nadgrids=conus,null \
    +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 44
-111 55
EOF
111d0'2.788"W   43d59'59.725"N 0.000
111dW   55dN 0.000

有关更多信息,请参阅有关的章节 其他变换网格 .

告诫

  • 格网重叠的位置(如圆锥体和 ntv1_can.dat 例如)无论是否合适,都将使用为点找到的第一个。所以,比如说, +nadgrids=ntv1_can.dat ,Conus将导致加拿大的数据被用于美国北部的一些地区,即使Conus数据是被批准用于该地区的数据。仔细选择文件和文件顺序是必要的。在某些情况下,可能需要将边界跨越数据集预分割为加拿大点和美国点,以便对其进行适当的格网转换

  • 通过将PROJïu DEBUG环境变量设置为一个值,可以找到有关正在应用的网格偏移的其他详细信息。这将导致输出到stderr关于哪个网格用于移动点、加载的各种网格的边界等等