pyproj EPSG:4326 – How to Convert Degrees to Radians in pyproj

coordinatesepsglatitude longitudepyproj

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?

Best Answer

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.

Related Question