2

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?

  • 1
    Can't you just test for 4326 in an if and pass the values as-is? – Vince Jun 21 at 10:16
  • Sure I can but that would be special case handling where it shouldn't be a special case. EPSG:4326 is one SRS I know of that is geodetic, i. e. uses degrees as unit instead of meters. There might (and will!) be others (e. g. the ancient one which just uses Paris for the longitude 0° and is otherwise the same), and I am reluctant to open my own list of SRSs now which should be handled specially. – Alfe Jun 21 at 12:37
  • This is perilously close to a bug report. Even if all the other GCS values are converted to radians, there isn't anything we can do about it – Vince Jun 21 at 13:33
  • I don't think inverse does what you think it does. It seems to just determine whether the transformation is from x/y→lon/lat or from lon/lat→x/y. Lon/lat will be degrees, and x/y will be metres (if using a projection with metres as a unit). Do a search for inverse here: pyproj4.github.io/pyproj/stable/_modules/pyproj/… – alphabetasoup Jun 25 at 4:09
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.

  • I think I see a road to my solution in your answer, it just isn't completely clear yet. Please keep in mind that my input is x/y in some given projection, so it can be degrees or meters or whatever. My output is supposed to be degrees always. Please adjust your answer accordingly to this case. – Alfe Jun 24 at 8:34
  • I now use Transformer.from_crs(given_srs, CRS('epsg:4326')).transform(x, y), and this seems to produce degree values for both geodetic SRSs (like EPSG:4326) and for meters-based ones (like EPSG:3857). Thanks for pointing me to this (hopefully) correct use of the pyproj library! If you adjust your answer to include this, I will accept it :-) – Alfe Jun 24 at 8:45
  • Glad to hear that it helps! I just updated it to better match your use case. – snowman2 Jun 25 at 2:10
1

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.

  • I'm missing a clear connection to my answer. Sure, I can guess on what you advise me to do now in my case (convert degrees/radians myself whenever I have a geographic projection maybe?), but a clear message would be helpful. – Alfe Jun 24 at 8:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.