0

I have about 20 million sets coordinates from the Philadelphia, PA, USA region. For each set of coordinates I would like the FIPS code of the corresponding census tract.

I looked at this thread: Free api to reverse geocode latitude, longitude to census tract? However, I think any http/API solution would be too slow. Also, I don't think QGIS could handle so many points.

Ideally I'd like to use Python for my solution. Maybe there is a function that can import polygons from a shp file and identify which contains a given latitude/longitude pair?

  • What spatial Python libraries have you tested (or at least looked into) so far? – PolyGeo Dec 8 '15 at 3:49
  • After some internet research/minimal ipython testing I am considering pyshp, pyproj and shapely. I'd also like to check out the tools mentioned by GrantD71. Any feedback on those? – Selah Dec 8 '15 at 16:08
1

If you want this to be fast and you have a large number of lat/longs to check I highly recommend using postgresql/PostGIS. However it is possible to do this in python, provided you have OGR/GDAL packages installed properly.

Here is an outline of how you might do this purely in python, assuming you have OGR/GDAL packages properly installed (these installations can be non-trivial, so beware)

from osgeo import ogr

#Assumes your points and shapefile are already in the same datum/projection
shapefile_name = "census_tracts.shp"

#This version takes a long_lat_list of the form below and a shapefile name
def getCensusTracts(long_lat_list, shapefile_name):
    driver = ogr.GetDriverByName("ESRI Shapefile")
    dataSource = driver.Open(shapefile_name, 0)
    layer = dataSource.GetLayer()
    results_dict = {}
    i = 0
    for feature in layer:
        geom = feature.GetGeometryRef()
        i += 1
        for pt in long_lat_list:
            gid = pt[0]
            lon = pt[1]
            lat = pt[2]
            point = ogr.Geometry(ogr.wkbPoint)
            point.AddPoint(lon, lat)
            if point.Within(geom) == True:
                feat_id = feature.GetField("fips")
                if gid in results_dict and feat_id not in results_dict[gid]:
                    results_dict[gid].append(feat_id)
                else:
                    results_dict[gid] = [feat_id]
    for pt in long_lat_list:
        gid = pt[0]
        lon = pt[1]
        lat = pt[2]
        if gid not in results_dict:
            results_dict[gid] = ['NA']
    return results_dict


#Where elements are [id, long, lat]
long_lat_list = [[1, -87.5, 35.5],[2, -78.446, 41.353]]

results_dict = getCensusTracts(long_lat_list, shapefile_name)
#results_dict returns a dictionary where {'id: list_of_fips_codes}
print results_dict

Edit: Given Selah's comments I discovered that you cannot iterate over a layer multiple times in OGR. I rewrote the solution so that the census shapefile only need be iterated over once.

  • Wierd... I get proper results for the first one or two sets of coordinates. After that I get only NA's. When I debug I get "Process finished with exit code 139" It seems like osgeo is crashing? I'm using version 1.10.1, will try updating versions. – Selah Dec 8 '15 at 20:32
  • Hmmm weird, that is a memory error. You might be running out of memory on your machine? I would check to see if you are going into swap (loading all census tracts into memory may not be an ideal solution for all machines, hence my postgresql suggestion) – GrantD71 Dec 8 '15 at 20:38
  • I have 30 gig of memory on my machine and was running on only 1000 pairs, so I wonder if it could be something different? – Selah Dec 8 '15 at 20:41
  • It's hackish, but I noticed if I reassign my, driver, dataSource and layer variables in the getCensusTract function I DO get results for every input. Might be good enough for now! – Selah Dec 8 '15 at 20:41
  • We use MySQL at my office so I was shying away from a solution that required migrating to postgresql... but I will keep it in mind. – Selah Dec 8 '15 at 20:42
0

I was using the FCC's API for a bit but for 19k coordinates, the process took several hours. GrantD71's script seemed to work pretty quick on a few hundred but I had issues getting my set done in a timely manner. Geopandas spatial join is the quickest way I've found so far. The TIGER shapefiles are on the Census's website. Import them and your coordinates like the below (you may have to project into different CRS - coordinate reference system).

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point
import os

os.environ['PROJ_LIB'] = r'C:\Users\root\Anaconda3\Library\share\proj'

census_tracts = gpd.read_file(r'C:\Users\root\Downloads\tl_2018_06_bg\tl_2018_06_bg.shp')

points_df = pd.read_csv(r'''C:\Users\root\Desktop\Desktop\airbnb_LAarea.csv''', index_col=0)

geometry = [Point(xy) for xy in zip(points_df.lng, points_df.lat)]
crs = {'init' :'epsg:4326'}
gdf = gpd.GeoDataFrame(points_df, crs=crs, geometry=geometry)

merged_file = gpd.sjoin(gdf, census_tracts, how='left', op='within')

merged_df = pd.DataFrame(merged_file)

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.