pyproj EPSG:4326 converts degrees to radians

Question

Case

My input is a point given as x, y and an SRS name which specifies how to interpret the numbers. My goal is to convert these into lon, lat in degrees.

I convert this by using this piece of code:

projection = pyproj.Proj(srs_name)
lon, lat = projection(x, y, inverse=True)

A typical input is coordinates in meters and the name of an SRS which has meters as unit. The code above works fine in this case.

Problem

Sometimes the input is the SRS name "EPSG:4326", and the input coordinates then are given in degrees. From what I understand about EPSG:4326 this is correct, too. Since the input already is in degrees, the conversion done by pyproj should not change anything, i. e. projection(x, y, inverse=True) should return a lon, lat with lon == x and lat == y.

Unfortunately this is not the case.

I figured by now that the projection = pyproj.Proj('EPSG:4326') converts the input from degrees to radians. So if I use the inverse=True variant, it converts radians to degrees. Since my input is not given in radians, the result will be nonsense (in degrees).

What I Tried

I played around with the kwargs radians=False and preserve_units=True. I tried giving it to the pyproj.Proj() and to the call of the conversion (projection(x, y, invers=True, radians=False)). None of the combinations I tried had any influence. The inverse=True conversion always assumes radians as input.

I also tried figuring out if there is some kind of global switch in pyproj to make it use degrees instead of radians. I didn't find any.

I'm using pyproj 2.2.0. From the git log of the project (as in Github) I understand that there have been some changes made in the recent past concerning the kwargs radians. Version 2.2.0 could be older than these, and the changes could be a fix for my case.

I also found on https://epsg.io/4326 that the unit of EPSG:4326 is specified as UNIT["degree",0.0174532925199433, ...] and this number is the factor between radians and degrees. This now looks like the unit of this EPSG in fact is radians. Should this be the case, then I wonder why my input is in degrees instead of radians if the SRS is EPSG:4326. Maybe the client (providing the input) does it wrong?

Question

How do I have to change my use of pyproj in order to not run into the problem anymore?

Answer 1

From the pyproj FAQ page, it recommends using the Transformer class for this type of operation as it handles datum shifts. I would recommend reading http://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-longitude-to-projection-converter as it will help clarify if you should use the geodetic CRS or EPSG:4326 in your use case.

So, for your use case, you would likely want to use:

>>> from pyproj import CRS, Transformer
>>> from_crs = CRS("epsg:4326") # or whatever you want
>>> transformer = Transformer.from_crs(from_crs, "epsg:4326")
>>> transformer.transform(-87.93514385942143, 42.000889116)
(-87.93514385942143, 42.000889116)
>>> from_crs = CRS("epsg:3857")
>>> transformer = Transformer.from_crs(from_crs, "epsg:4326")
>>> transformer.transform(4675517.589257867, -25615907.92738164)
(-87.93514385942143, 42.000889116)

Also, if you want your input to always be in the x,y order:

>>> from pyproj import Transformer
>>> transformer = Transformer.from_crs("epsg:26917", "epsg:4326", always_xy=True)
>>> transformer.transform(571666.4475041276, 5539109.815175673)
(-80.0, 50.0)

This is useful as the axis order is not always consistent. See: https://proj.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent

Hopefully this is helpful and will get you headed in the right direction.

Answer 2

I had a hard time tracking this down but finally looked in the pyproj tests for CRS.

What you need to do is test whether the input CRS is geographic or not. Test is directly from pyproj python script (above link).

def test_is_geographic():
    assert CRS({"init": "EPSG:4326"}).is_geographic is True
    assert CRS({"init": "EPSG:3857"}).is_geographic is False

    wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
    assert wgs84_crs.is_geographic is True

    nad27_crs = CRS.from_string("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs")
    assert nad27_crs.is_geographic is True

    lcc_crs = CRS.from_string(
        "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 "
        "+units=m +lat_2=77 +lat_1=49 +lat_0=0"
    )
    assert lcc_crs.is_geographic is False

You also might want to look at using Transformer or double-check what your code is doing now. "Proj" in pyproj says that it's converting to or from the geographic CRS that's used in the Projected CRS. If the ProjCRS isn't based on WGS84, you're not getting back WGS84 latitude-longitude values.