从 CRS 创建转换器时出错:源椭球体和目标椭球体不属于同一个天体

Question

我正在尝试使用pyproj的CRS创建一个转换。我想将立体投影中的地图数据(荷兰(转换为纬度经度表示。我在地图的元数据中找到了必要的转换信息，但出现以下错误：

pyproj.exceptions.ProjError: Error creating Transformer from CRS.: (Internal Proj Error: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body)

我使用以下 python3 代码生成转换器：

import h5py
from pyproj import CRS, Transformer
with h5py.File(folder+filename, 'r') as f:
proj4 = str(list(f['geographic/map_projection'].attrs.items())[2][1])
proj4 = proj4[2:len(proj4)-1]
from_proj = CRS.from_proj4(proj4)
to_proj = CRS.from_proj4("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
print(from_proj)
print(to_proj)
transform = Transformer.from_crs(from_proj, to_proj)

print(from_proj)

输出：

+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378.137 +b=6356.752 +x_0=0 +y_0=0 +type=crs

print(to_proj)

输出：

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 +type=crs

transform = Transformer.from_crs(from_proj, to_proj)

生成错误：

Traceback (most recent call last):
File "load_random_h5.py", line 100, in <module>
load_random_h5()
File "load_random_h5.py", line 50, in load_random_h5
transform = Transformer.from_crs(from_proj, to_proj)
File "/Users/user/miniconda2/envs/gdal_test/lib/python3.7/site-packages/pyproj/transformer.py", line 323, in from_crs
area_of_interest=area_of_interest,
File "pyproj/_transformer.pyx", line 311, in pyproj._transformer._Transformer.from_crs
pyproj.exceptions.ProjError: Error creating Transformer from CRS.: (Internal Proj Error: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body)

问题可能源于这样一个事实，即地图数据表示空气中约1公里的降水(水落在大气中(，因此与地球半径不匹配，但这只是猜测。地图数据是荷兰两个雷达站的复合数据。但我认为元数据中的投影信息应该足以应用转换。

我试图用 EPSG 标准的几个投影替换from_proj投影，但它们都返回无意义的经度纬度坐标(在荷兰的情况下，感觉是纬度在 50 到 53 之间，经度在 3 到 8 之间(。将标志+ellps=WGS84 +towgs84=0,0,0连接到proj4可以消除错误，但会再次返回无意义的经度纬度坐标。

有没有人知道解决错误的方法或修复它的方法？

Answer 1

我相信投影中定义的a和b参数的单位不正确。当它们需要以米为单位时，它们似乎以公里为单位。

当查看+proj=latlon的椭球体参数时，您可以看到椭球体的星等是其他投影的~1000倍：

>>> from pyproj import CRS
>>> cc = CRS("+proj=longlat")
>>> cc.datum.ellipsoid.semi_major_metre
6378137.0
>>> cc.datum.ellipsoid.semi_minor_metre
6356752.314245179
>>> cp = CRS("+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378.137 +b=6356.752 +x_0=0 +y_0=0 +type=crs")
>>> cp.datum.ellipsoid.semi_major_metre
6378.137
>>> cp.datum.ellipsoid.semi_minor_metre
6356.752

基于此，将a和b参数乘以 1000 应该可以修复您的椭球体：

>>> cp = CRS("+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378137 +b=6356752 +x_0=0 +y_0=0 +type=crs")
>>> cp.datum.ellipsoid.semi_major_metre
6378137.0
>>> cp.datum.ellipsoid.semi_minor_metre
6356752.0

创建转换器时没有错误：

>>> from pyproj import Transformer
>>> trans = Transformer.from_crs("+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378137 +b=6356752 +x_0=0 +y_0=0 +type=crs", "+proj=latlon")
>>> trans
<Concatenated Operation Transformer: pipeline>
Description: Inverse of unknown + Ballpark geographic offset from unknown to unknown
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)