解决pyproj高程基准转换失败问题

解决pyproj高程基准转换失败问题

什么是高程基准转换？

以新西兰坐标系为例，其所采用的大地坐标系为NZGD2000（同我们熟知的WGS84、CGCS2000），高程起算面为GRS80椭球面，也就是图中的NZGD2000 Ellipsoid，这是一个数学上定义的曲面，以该曲面起算的高程为椭球高HAE。而我们工程中通常使用的高程是通过水准联测得到的高程，显然不是以椭球面作为高程起算面的，而是理论上以大地水准面作为高程起算面，但大地水准面又是不可求得的。因此各地均采用通过长期验潮站确定的平均海平面作为当地的高程起算面。这就出现一个问题，各地确定的平均海平面是不一样的，也就无法进行全球统一。因此在工程应用中，我们通常需要进行高程基准的转换，例如将椭球高转换为水准高（即以平均海平面作为起点的高程，Mean Sea Level，MSL）。也就是将地面点到椭球面NZGD2000 Ellipsoid的高差h_NZGD2000 转换为地面点到平均海平面NZGeoid06的高差H_NZVD2016 。

利用Python进行高程基准转换

例如，将WGS84椭球高转换为NZVD2016高程，可通过以下代码实现：

from pyproj import Transformer
transformer_latlon = Transformer.from_crs(  # 高程转换"epsg:4326", "epsg:7839")  # WGS84->NZVD2016
lat, lon, hei = transformer_latlon.transform( # lat，lon，hei为转换后坐标lat0, lon0, hei0)  # lat0，lon0，hei0为转换前坐标

准换完成后，可以看一下高程值是否发生了变化，如果有变化说明转换成功后面内容就不用看了。
如果没有，那就是转换失败，需要后面继续解决问题了。

解决高程基准转换失败问题

前面说过，进行高程基准准换时需要当地的平均海平面的，因此转换失败很大可能就是因为这个当地平均海平面数据缺失导致的。那么，如何确定真的是数据缺失呢？方法如下：

PS E:\Code\VSCode> python # 命令行打开python
Python 3.8.0 (tags/v3.8.0:fa919fd, Oct 14 2019, 19:37:50) [MSC v.1916 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.>>> from pyproj.transformer import TransformerGroup # 导入 pyproj
>>> tg = trans_group = TransformerGroup(4326, 7839)
C:\Users\Dell\AppData\Local\Programs\Python\Python38\lib\site-packages\pyproj\transformer.py:197: UserWarning: Best transformation is not available due to missing Grid(short_name=nz_linz_nzgeoid2016.tif, full_name=, package_name=, url=https://cdn.proj.org/nz_linz_nzgeoid2016.tif, direct_download=True, open_license=True, available=False)super().__init__(
>>> tg
<TransformerGroup: best_available=False>
- transformers: 0
- unavailable_operations: 1

如上可见输出“Best transformation is not available due to missing Grid” ，表明缺少相关的水准栅格数据。
那么就需要单独针对这个基准进行下载了：
我们先来看一下缺少栅格的信息：

>>> tg.transformers[0].description
Traceback (most recent call last):File "<stdin>", line 1, in <module>
IndexError: list index out of range
>>> tg.unavailable_operations[0].name
'NZGD2000 to NZVD2016 height (2)'
>>> tg.unavailable_operations[0].grids[0].url
'https://cdn.proj.org/nz_linz_nzgeoid2016.tif'

可以看到是可以找到对应tif的下载路径的，这表明可以自动进行下载：

>>> tg.download_grids(verbose=True)
Downloading: https://cdn.proj.org/nz_linz_nzgeoid2016.tif

当找不到这个路径时，说明路径缺失了，需要进一步操作，具体参考：link

参考：

• https://pyproj4.github.io/pyproj/stable/transformation_grids.html