This is a Python question, involving geopandas and the libraries it calls, primarily pandas and shapely.
tl;dr:
I am trying to run a test using
if foo.intersects(bar): and am repeatedly getting one of two errors:
ValueError: The truth value of a Series is ambiguous. Use a.empty, a.bool(), a.item(), a.any() or a.all().AttributeError: 'GeoSeries' object has no attribute '_geom'
I believe the problem is that I am running this boolean intersects() test on a shapely geometry object and a GeoSeries. I cannot figure out how to access the geometry of a single item within a GeoSeries.
More details
I have a list of >100,000 establishments located all over the United States. I need to calculate a dyadic measure for pairs of establishments that are within a certain distance (about a quarter-mile) of each other.
I created two GeoDataFrames, one with the point list of establishments, the other a polygon list of circles centered on each establishment, with a radius of one quarter mile. Within each GeoDataFrame, I also know the county FIPS within which dot is located (or the circle is centered).
I wrote code that iterates over each circle and uses .within() to return a GeoDataFrame with just the establishments inside the circle. This works but is incredibly slow, because it has to test every establishment for every circle.
Given the small size of the circles, I realized I only really needed to check establishments in the county where the focal establishment was located, plus (to be certain) adjacent counties. (I have a shapefile with county geometries.) I put together a dictionary with each county's adjacent counties, and modified my code to check all establishments in that subset of the original list. This works, and is faster, but still not fast.
I then thought that I should first check which counties the focal circle intersects, and only do the tests on establishments in the relevant counties. This is where I am now stuck.
Below I provide code, including my first two approaches and the current one that isn't working:
import pandas
import geopandas
import os.path
from counties_adjacencies import county_adjacent
# county_adjacent is a dictionary containing, for each county FIPS
# code, the adjacent counties' FIPS codes. That is:
# county_adjacent = {
# 1001 : [1001,1021,1047,1051,1085,1101],
# 1003 : [1003,1025,1053,1097,1099,1129,12033],
# 1005 : [1005,1011,1045,1067,1109,1113,13061,13239,13259],
# ...
# 78020 : [78020,78030],
# 78030 : [78020,78030]
# }
# I have two ESRI shapefiles, 'ests_1971.shp' and
# 'ests_1971_buffer.shp'. The first is a point file with establishment
# locations, the second is a polygon file with 1/4-mile-radius circles
# centered on each establishment. I also have the TIGER line shape
# file with county geometries.
circles = geopandas.read_file(
os.path.relpath("../data/ests_1971_buffer.shp"))
establishments = geopandas.read_file(
os.path.relpath("../data/ests_1971.shp"))
counties = geopandas.read_file(
os.path.relpath("../canonical/tl_2017_us_county.shp"))
# I add a county FIPS code in the county shape file that has the
# state FIPS but no leading zeros, to match the county FIPS
# codes in the other shape files.
counties['county'] = ((counties['STATEFP'] +
counties['COUNTYFP']).str.lstrip('0')).astype(int)
# Create a dataframe into which I'll put the list of relevant dyads.
dyads = pandas.DataFrame()
# FIRST APPROACH: For each circle, test all establishments to see
# whether they lie inside the circle. This works but is appallingly
# slow, because it has to look at many, many establishments that
# aren't really candidates for lying inside.
for circle, estloc, estid, estcounty in zip(circles.geometry,
establishments.geometry,
establishments.estid,
establishments.county):
subests = establishments[establishments.within(circle)].copy()
# DO SOME STUFF USING ESTLOC AND ESTID HERE ...
# ...
dyads = pandas.concat([dyads, subests])
# SECOND APPROACH: For each circle, test all establishments within the
# circle's county and its adjacent counties to see whether they lie
# inside the circle. This also works but is still pretty slow, because
# it tests the establishments in counties that may lie outside the
# circle radius.
for circle, estloc, estid, estcounty in zip(circles.geometry,
establishments.geometry,
establishments.estid,
establishments.couty):
for county in county_adjacent[estcounty]:
county_ests = establishments[establishments['county'] == county]
subests = county_ests[county_ests.within(circle)].copy()
# DO SOME STUFF USING ESTLOC AND ESTID HERE ...
# ...
dyads = pandas.concat([dyads, subests])
# THIRD APPROACH: For each circle, test which counties it intersects,
# then test establishments within the intersecting counties to see
# whether they lie inside the circle. In the vast majority of cases,
# this will only require testing establishments within a single
# county, so should run much faster. But I can't get the intersection
# test to work.
for circle, estloc, estid, estcounty in zip(circles.geometry,
establishments.geometry,
establishments.estid,
establishments.county):
for county in county_adjacent[estcounty]:
the_county = counties[counties['county'] == county].geometry
if the_county.intersects(circle): # THIS IS WHERE THINGS BREAK
county_ests = establishments[establishments['county'] == county]
subests = county_ests[county_ests.within(circle)].copy()
# DO SOME STUFF USING ESTLOC AND ESTID HERE ...
# ...
dyads = pandas.concat([dyads, subests])
else:
continue
Note that I'm trying to run .intersects() on circle, which is a shapely geometry object, and the_county, which is a GeoDataFrame (with one row, but is a GeoSeries nonetheless). I suspect there is a way to get this intersection test to work, but I'm running out of ideas.
Can someone point out what I'm doing wrong here?
