1

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.

2

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).

  • You're right, thanks. I've switched my latitude and longitude The incorrect projection was an artifact of testing that I have also corrected. The issue still persists. The problem is that this projection is scaling down my latitude and longitude to around the origin which is nowhere near the point I'm trying to plot. Also, an interesting thing I found is that if I do not project the point using pyproj it ends up basically where I want it plotted. Why would projecting the point using the same projection as the map cause such distortion? – Allen More Aug 28 '17 at 15:19
  • This is no longer the same question. – BradHards Aug 29 '17 at 3:47
1

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)

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.