Unexpected Conversion Result using Pyproj

Question

I am trying to use Pyproj to convert a location in Boston from State Plane to Latitude and Longitude. However, the result did not turn out as I expected. Here is my snippet:
from pyproj import Proj, transform

inProj = Proj(init='epsg:2229', preserve_units=True)
outProj = Proj(init='epsg:4326')
x1,y1 = 783330.132997244, 2963345.94232371
LONGITUDE,LATITUDE = transform(inProj,outProj,x1,y1)
print(LATITUDE,LONGITUDE)

The result is "35.5856776921177 -137.5552510169851" which is way too far away from the correct one calculated via https://www.earthpoint.us/StatePlane.aspx?fbclid=IwAR284hbFbNerEsCcvKtJPCG20eOAk5ft8oEdAi3KGRUQLqEXsLtmRHQb40s.
Note: my data is given in US Survey Feet and Zone Number should be 2001-Massachusetts Mainland.

snowman2 · Answer

There are two things to watch out for:

The CRS you specified is meant for California:

>>> from pyproj import CRS
>>> CRS("epsg:2229")
<Projected CRS: EPSG:2229>
Name: NAD83 / California zone 5 (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - California - counties Kern; Los Angeles; San Bernardino; San Luis Obispo; Santa Barbara; Ventura.
- bounds: (-121.42, 32.76, -114.12, 35.81)
Coordinate Operation:
- name: SPCS83 California zone 5 (US Survey feet)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

I think you want something more like:

>>> CRS("EPSG:2249")
<Projected CRS: EPSG:2249>
Name: NAD83 / Massachusetts Mainland (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - Massachusetts onshore - counties of Barnstable; Berkshire; Bristol; Essex; Franklin; Hampden; Hampshire; Middlesex; Norfolk; Plymouth; Suffolk; Worcester.
- bounds: (-73.5, 41.46, -69.86, 42.89)
Coordinate Operation:
- name: SPCS83 Massachusetts Mainland zone (US Survey feet)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

Be careful of axis order. The order is actually latitude, longitude for EPSG:4326 (only if not using init=). You can make it longitude, latitude with always_xy=True

With both of those combined with https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1, your solution would likely come from:

>>> from pyproj import Transformer
>>> transformer = Transformer.from_crs("epsg:2249", "epsg:4326", always_xy=True)
>>> transformer.transform(783330.132997244, 2963345.94232371)
(-71.02936000125914, 42.37861000033115)

Unexpected Conversion Result using Pyproj

One Answer

Add your own answers!

Ask a Question