Convert coordinates to census tract for large dataset

Question

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?

Answer 1

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)

Answer 2

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.