Something going wrong trying to transform longitude, latitude to NAD 83 using pyproj

Question

I'd like to convert a latitude longitude pair to NAD 83. I've read NAD 83 is equivalent to epsg:4269 and the shapefile I am working with is a US Census TIGER shapefile, all of which are supposed to be in NAD 83 and geopandas reports my GeoDataFrame as being in epsg 4269. All of this makes sense. The issue I'm running into occurs when I want to plot a point described as a lat/long pair onto the shapefile. My understanding is that, in order to do this, I need to transform the point using the relevant map projection (NAD 83/epsg:4269), however, when I try to do this using pyproj, the point returned is way off:

import geopandas as gp
import matplotlib.pyplot as plt
from shapely.geometry import Point
import pyproj

gdf = gp.read_file('./state_shapefiles/WI/tl_2016_55_bg.shp')
projection = pyproj.Proj(init='epsg:4269')

lat = 43.087806
lon = -89.430121


# this returns tuple(-1.5608500619069237, 0.7520240821605678) which is 
# nowhere near the point I am trying to indicate on my map (Madison
proj_coords = projection(lon, lat)

point = Point(proj_coords[1], proj_coords[0]).buffer(1.0)

layer = gp.GeoDataFrame({'geometry': [point]})

layer_plot = layer.plot(color='red')

gdf.plot(ax=layer_plot)
plt.show()

print('end')

Am I misinterpreting something or have implemented this wrong?

Edit: Fixed latitude and longitude mix up, issue persists.

Answer 1

Probably that lon = 43.087806 lat = -89.430121 is not what you wanted (latitude -89 degrees is near the south pole).

Also, it isn't clear why you are initialising to EPSG:3857 - you need to use the real coordinate basis (e.g. EPSG:4269).

Answer 2

I would recommend using the pyproj.Transformer class as it also takes into account datum shifts (See: https://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-longitude-to-projection-converter):

>>> from pyproj import Transformer
>>> transformer = Transformer.from_crs("epsg:4326", "epsg:4269", always_xy=True)
>>> transformer.transform(-89.430121, 43.087806)
(-89.430121, 43.087806)

Although the point does not appear to have shifted much in this case.

>>> from pyproj import CRS
>>> 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

>>> crs_4269 = CRS("epsg:4269")
>>> crs_4326 = CRS("epsg:4326")
>>> crs_4326.datum == crs_4269.datum
False
>>> crs_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

>>> crs_4326
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

For reference:

$ python -c "import pyproj; pyproj.show_versions()"

System:
    python: 3.7.3 | packaged by conda-forge | (default, Mar 27 2019, 23:01:00)  [GCC 7.3.0]
executable: ~/miniconda3/envs/pyproj/bin/python
   machine: Linux-4.15.0-51-generic-x86_64-with-debian-buster-sid

PROJ:
      PROJ: 6.1.1
  data dir: ~/scripts/pyproj/pyproj/proj_dir/share/proj

Python deps:
    pyproj: 2.2.1
       pip: 18.0
setuptools: 41.0.1
    Cython: 0.29.10
     aenum: (2, 1, 2)