pyproj in Python returning incorrect coordinates for NAD_1927_StatePlane_Wisconsin_South_FIPS_4803

Question

I'm new to using geospatial data. I'm having trouble with pyproj getting back the correct longitude and latitude from a file with a different cooridinate system. I'm working with a file that has the Milwaukee city limits that is in State Plane coordinates NAD_1927_StatePlane_Wisconsin_South_FIPS_4803 and I want the WGS84 coordinates, which should be just the ordinary longitude and latitude. I want to convert this because it would be easier work with tools like folium that I'm using.

https://data.milwaukee.gov/dataset/corporate-boundary The first coordinates in the shp file are [2516440.160073474, 440481.9301029742]

It has the prj file:

PROJCS["NAD_1927_StatePlane_Wisconsin_South_FIPS_4803",
       GEOGCS["GCS_North_American_1927",
              DATUM["D_North_American_1927",
                    SPHEROID["Clarke_1866",6378206.4,294.9786982]],
              PRIMEM["Greenwich",0.0],
              UNIT["Degree",0.017453292519943295]],
       PROJECTION["Lambert_Conformal_Conic"],
       PARAMETER["False_Easting",2000000.0],
       PARAMETER["False_Northing",0.0],
       PARAMETER["Central_Meridian",-90.0],
       PARAMETER["Standard_Parallel_1",42.73333333333333],
       PARAMETER["Standard_Parallel_2",44.06666666666667],
       PARAMETER["Latitude_Of_Origin",42.0],
       UNIT["Foot_US",0.30480060960121924]]

According to https://spatialreference.org this would have the code EPSG6948, but pyproj doesn't recognize that so I have build a custom string. After looking at documentation and examples I came up with the code:

import pyproj
from pyproj import Transformer
cust = "+proj=lcc +ellps=clrk66 +datum=NAD27 +lat_1=42.73333333333333 +lat_2=44.06666666666667 +lat_0=42.0 +lon_0=-90.0 +x_0=2000000.0 +y_0=0.0 +units=us-ft +no_defs"
t = Transformer.from_crs(cust, "WGS84")
x1,y1 = [2516440.160073474, 440481.9301029742]
lat, lon = t.transform(x1, y1)
print(lat, ", ", lon)

According to the site https://www.earthpoint.us/StatePlane.aspx this code should return the coordinates: Degrees Lat Long 43.1900328, -087.9451397

But the code prints out: 42.206953018562416 , -105.00939853775037 which is somewhere in Wyoming not Milwaukee Wisconsin. I've check over this a few times and can't see what I'm doing wrong. If anyone can help I would much appreciate it.

snowman2 snowman2 · Accepted Answer · 2021-03-26T00:47:26

You can use the WKT string to create the transformer:

>>> from pyproj import Transformer
>>> wkt = """PROJCS["NAD_1927_StatePlane_Wisconsin_South_FIPS_4803",
...        GEOGCS["GCS_North_American_1927",
...               DATUM["D_North_American_1927",
...                     SPHEROID["Clarke_1866",6378206.4,294.9786982]],
...               PRIMEM["Greenwich",0.0],
...               UNIT["Degree",0.017453292519943295]],
...        PROJECTION["Lambert_Conformal_Conic"],
...        PARAMETER["False_Easting",2000000.0],
...        PARAMETER["False_Northing",0.0],
...        PARAMETER["Central_Meridian",-90.0],
...        PARAMETER["Standard_Parallel_1",42.73333333333333],
...        PARAMETER["Standard_Parallel_2",44.06666666666667],
...        PARAMETER["Latitude_Of_Origin",42.0],
...        UNIT["Foot_US",0.30480060960121924]]"""
>>> transformer = Transformer.from_crs(wkt, "WGS84")
>>> transformer.transform(2516440.160073474, 440481.9301029742)
(43.19212933447268, -88.06334977811528)

Side note, I think this is the EPSG code you want:

<Projected CRS: EPSG:32054>
Name: NAD27 / Wisconsin South
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - Wisconsin - counties of Adams; Calumet; Columbia; Crawford; Dane; Dodge; Fond Du Lac; Grant; Green; Green Lake; Iowa; Jefferson; Juneau; Kenosha; La Crosse; Lafayette; Manitowoc; Marquette; Milwaukee; Monroe; Ozaukee; Racine; Richland; Rock; Sauk; Sheboygan; Vernon; Walworth; Washington; Waukesha; Waushara; Winnebago.
- bounds: (-91.43, 42.48, -86.95, 44.33)
Coordinate Operation:
- name: Wisconsin CS27 South zone
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1927
- Ellipsoid: Clarke 1866
- Prime Meridian: Greenwich

>>> transformer = Transformer.from_crs("EPSG:32054", "WGS84")
>>> transformer.transform(2516440.160073474, 440481.9301029742)
(43.19212933447268, -88.06334977811528)

pyproj in Python returning incorrect coordinates for NAD_1927_StatePlane_Wisconsin_South_FIPS_4803

1 Answers