python pyproj transform yielding different results for the same input parameters

Question

I'm running a simple python script to transform longlat Nad83 to Ky State Plane South Nad83 (epsg 4269 to 6474). When I set the input Projection to:

inProj = pyproj.Proj(proj='longlat', datum ='NAD83')

I get one answer (which appears to be correct). When I set the input projection :

inProj = pyproj.Proj('epsg:4269')

I get a different answer. I print out inProj both times and both appear to be the same.

Run 1:

import pyproj
import numpy as np

# nad 83 latt and longs      
inProj        = pyproj.Proj('epsg:4269')      

#  KY State Nad83 Meter (2011)
outProj  = pyproj.Proj("epsg:6474") #  KY State South Nad83 Meter (2011)

long = np.array([-83.88259259,-83.88259259,-83.55083333,-83.55083333])
latt = np.array([ 36.49333333, 36.79203704, 36.79203704, 36.49333333])

x,y = pyproj.transform(inProj, outProj ,long ,latt)

print('\ninProj = ',inProj,'\n\noutProj = ',outProj,'\n')

for i in range(len(long)):
    print('long = ',long[i],'latt = ',latt[i],
          'x = ',x[i],'y = ',y[i])

Output:

inProj =  Proj('+proj=longlat +datum=NAD83 +no_defs', preserve_units=True) 

outProj =  Proj('+proj=lcc +lat_0=36.3333333333333 +lon_0=-85.75 +lat_1=37.9333333333333 +lat_2=36.7333333333333 +x_0=500000 +y_0=500000 +ellps=GRS80 +units=m +no_defs', preserve_units=True) 

long =  -83.88259259 latt =  36.49333333 x =  72927331.14729036 y =  -11599224.291812904
long =  -83.88259259 latt =  36.79203704 x =  72992044.27873956 y =  -11370127.68625123
long =  -83.55083333 latt =  36.79203704 x =  70702668.24198839 y =  -10727351.707669286
long =  -83.55083333 latt =  36.49333333 x =  70639998.82021627 y =  -10949213.198204307

Run2:

import pyproj
import numpy as np

# nad 83 latt and longs   
inProj        = pyproj.Proj(proj='longlat', datum ='NAD83')     


#  KY State Nad83 Meter (2011)
outProj  = pyproj.Proj("epsg:6474") #  KY State South Nad83 Meter (2011)

long = np.array([-83.88259259,-83.88259259,-83.55083333,-83.55083333])
latt = np.array([ 36.49333333, 36.79203704, 36.79203704, 36.49333333])

x,y = pyproj.transform(inProj, outProj ,long ,latt)

print('\ninProj = ',inProj,'\n\noutProj = ',outProj,'\n')

for i in range(len(long)):
    print('long = ',long[i],'latt = ',latt[i],
          'x = ',x[i],'y = ',y[i]

Output:

inProj =  Proj('+proj=longlat +datum=NAD83 +no_defs', preserve_units=True) 

outProj =  Proj('+proj=lcc +lat_0=36.3333333333333 +lon_0=-85.75 +lat_1=37.9333333333333 +lat_2=36.7333333333333 +x_0=500000 +y_0=500000 +ellps=GRS80 +units=m +no_defs', preserve_units=True) 

long =  -83.88259259 latt =  36.49333333 x =  667315.2681430576 y =  519409.6151513982
long =  -83.88259259 latt =  36.79203704 x =  666660.1039730347 y =  552551.2004212044
long =  -83.55083333 latt =  36.79203704 x =  696263.6061808041 y =  553188.4213983886
long =  -83.55083333 latt =  36.49333333 x =  697035.1458569702 y =  520049.34113333403

Pyproj appears to be treating the input parameters as the same but the results are much different.

snowman2 snowman2 · Accepted Answer · 2020-03-02T13:42:09

You are running into axis order changes introduced with PROJ 6:

Traditionally, the order has been east(x) north(y) as shown in the PROJ param version:

>>> from pyproj import CRS
>>> CRS("+proj=longlat +datum=NAD83")
<Geographic 2D CRS: +proj=longlat +datum=NAD83 +type=crs>
Name: unknown
Axis Info [ellipsoidal]:
- lon[east]: Longitude (degree)
- lat[north]: Latitude (degree)
Area of Use:
- undefined
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

But, in the EPSG version, it is North(y), East(x) order:

>>> CRS("EPSG:4269")
<Geographic 2D CRS: EPSG:4269>
Name: NAD83
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: North America - NAD83
- bounds: (167.65, 14.92, -47.74, 86.46)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

So, if you want it to give the correct output, either:

Switch the order of input:

>>> transformer = Transformer.from_crs("EPSG:4269", "epsg:6474")
>>> long, latt = -83.88259259, 36.49333333
>>> transformer.transform(latt, long)
(667315.2681430576, 519409.6151513982)

Use always_xy:

>>> transformer_xy = Transformer.from_crs("EPSG:4269", "epsg:6474", always_xy=True)
>>> transformer_xy.transform(long, latt)
(667315.2681430576, 519409.6151513982)

I would also recommend reading: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6

python pyproj transform yielding different results for the same input parameters

1 Answers