从CRS创建变压器时出错：源椭球和目标椭球不属于同一天体

Question

我正在尝试使用pyproj的CRS创建转换。我想将map data中的stereographic projection（荷兰）转换为纬度经度表示。我在地图的元数据中找到了必要的转换信息，但出现以下错误：

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)

[该问题可能源于以下事实：map data代表空气中约1公里处的降水（降落到大气中的水），因此与地球半径不匹配，但这只是推测。 map data由荷兰的两个雷达站组成。但是我认为元数据中的投影信息应该足以进行转换。

[我试图用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)