I am using Census Shapefiles from NHGIS. I am trying to convert points from it to easy to read latitude and longitude - but am having some trouble. The sample data I'm using for now is data from Alabama.

import geopandas
import random
from shapely.geometry import *
import pyproj
poly  = geopandas.GeoDataFrame.from_file('NewShapes/nhgis0103_shapefile_tl2010_010_block_2010/AL_block_2010.shp')
print poly.crs
p = pyproj.Proj("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
for index, row in poly.iterrows():
    print LY,BY,LX,BX
    lon, lat = p(BX, BY, inverse=True)
    print lon, lat

The output is:

{u'lon_0': -96, u'datum': u'NAD83', u'y_0': 0, u'no_defs': True, u'proj': u'aea', u'x_0': 0, u'units': u'm', u'lat_2': 45.5, u'lat_1': 29.5, u'lat_0': 37.5}
-515926.764688 -511955.403286 856239.343387 860999.442418
-88.180859033 17.9110274036

The latitude appears to be correct - but the longitude is way off. This is an example but I have tried other state shapefiles and other blocks within the data and getting similarly weird results.

Does anyone have an idea of what is going on here?

The bounds (-515926.764688 -511955.403286) do not appear to match any coordinate system I'm aware of.

  • The projection you defined (with p = pyproj.Proj(..)) does not match the crs of the dataframe you printed (eg lat_0 is different). You know there is a to_crs method on a GeoDataFrame with which you can also convert it to lat/lon? (poly.to_crs(epgs=4326).bounds) – joris Feb 7 '18 at 22:15
  • Thank you very much - I am still new to this but had tried for a couple hours with the manual solution and clearly made an error. This has been useful! – user7299543 Feb 7 '18 at 22:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.