CODEWITHSUNDEEP
python-3.xmap-projectionsproj
Christiaan Lamers
Christiaan Lamers

Reputation: 81

Error creating Transformer from CRS: Source and target ellipsoid do not belong to the same celestial body

I am trying to create a transformation using pyproj's CRS. I want to transform map data (of the Netherlands) in stereographic projection to a latitude longitude representation. I have found the necessary transformation info in the meta-data of the map, but I get the following error:

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)

I use the following python3 code to generate the transformator:

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)

outputs:

+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)

outputs:

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

transform = Transformer.from_crs(from_proj, to_proj)

generates an error:

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)

The problem might stem from the fact that the map data represents precipitation (water falling in the atmosphere) about 1 km in the air, thereby not matching the radius of the earth, but this is just speculation. The map data is a composite of two radar stations in the Netherlands. But I assume the projection information in the meta-data should be enough to apply a transform.

I have tried to replace the from_proj projection with several projections from the EPSG standard, but they all return nonsensical longitude latitude coordinates (sensical would be a latitude between 50 and 53, and a longitude between 3 and 8, in the case of the Netherlands). Concatinating the flags +ellps=WGS84 +towgs84=0,0,0 to proj4 negates the error, but again returns nonsensical longitude latitude coordinates.

Does anybody know a way around the error, or a way to fix it?

Upvotes: 1

Views: 3432

Answers (1)

snowman2
snowman2

Reputation: 711

I believe that the a and b parameters defined in your projection are in the incorrect units. They appear to be in kilometers when they need to be in meters.

When looking at the ellipsoid parameters of the +proj=latlon, you can see the magnitude of the ellipsoid is ~1000x that of the other projection:

>>> 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

Based on this, multiplying the a and b parameters by 1000 should fix your ellipsoid:

>>> 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

And it does not have an error when creating the transformer:

>>> 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)

Upvotes: 2

Related Questions