# Performing Spatial Join / match Points from dataframe to polygons using Python?

I have a dataframe with coordinates and other attributes, and a shp file (the whole package with shx and dbf as well) with many polygons of neighborhoods. I need to match each point to which polygon it belongs to, and it would be for quite a large dataset.

## What I did find and why it didn't fit

I've found some suggested solutions online for it, and understood the best ways should be by using Fiona and Shapely, or GeoPandas, but:

1. Existing solutions (that I found) were all assuming I'm working with an shp file for the points as well, but I'm working with a dataframe. (I may be able to find out how to convert it to an shp file but I assume that will be redundant and not the most efficient).

2. Most solutions assumed the shp file has only one attribute, while mine has many attributes; or aimed at checking if points are in 1 polygon or not. I assume there must be a better way to assign points to the corresponding polygon than looping through this.

## Solutions I tried

I tried suggestions from the following pages:

### Attempt 1: Fiona + Shapely

``````import fiona
from shapely.geometry import Point, shape

NY_nbr_shpfile = 'taxi_zones.shp'

def coor_to_nbr(longit, lat, shape_file):
mypoint = Point(longit, lat)
with fiona.open(NY_nbr_shpfile) as shp:
polygons = [poly for poly in shp]
poly_idx = [i for i, poly in enumerate(polygons)
if mypoint.within(shape(poly['geometry']))]
if poly_idx: print poly_idx
if not poly_idx:
return None
else:
# Take first polygon that overlaps since may overlap with several if on border
match = polygons[poly_idx[0]]
return match['properties']

### TRYING TO MATCH one-by-one, that didn't work:
print coor_to_nbr(-73.9868805930018, 40.7697167683218, NY_nbr_shpfile)
print coor_to_nbr( 40.722249, -73.997673, NY_nbr_shpfile) . ### just checking it's not opposite
print coor_to_nbr( -73.977673,40.722249, NY_nbr_shpfile)`
``````

This printed 'None' to everything I tried, but I know these points are inside one of the polygons.

### Attempt 2: GeoPandas1

This is taken from the following url; I left the original points there just to see if it would even work, and it seems to not work for me - yielding at ValueError: need at least one array to concatenate.

``````#### https://gis.stackexchange.com/questions/190903/assign-a-point-to-polygon-using-pandas-and-shapely ####
ny_nbr_shpfile = '/Users/tomer/Dropbox (Via)/Mapping data/New York/TLC Taxi Zones/taxi_zones (1)/taxi_zones.shp'

import pandas
import geopandas
import geopandas.tools
from shapely.geometry import Point

#### POLYGONS
polys = geopandas.GeoDataFrame.from_file(ny_nbr_shpfile)

#### POINTS
# Create a DataFrame with some cities, including their location
raw_data = [
("London", 51.5, -0.1),
("Paris", 48.9, 2.4),
("San Francisco", 37.8, -122.4),
]
points = pandas.DataFrame(raw_data, columns=["name", "latitude", "longitude"])
###print points
# Create the geometry column from the coordinates
# Remember that longitude is east-west (i.e. X) and latitude is north-south (i.e. Y)
points["geometry"] = points.apply(lambda row: Point(row["longitude"], row["latitude"]), axis=1)
del(points["latitude"], points["longitude"])
###print points
# Convert to a GeoDataFrame
points = geopandas.GeoDataFrame(points, geometry="geometry")

# Declare the coordinate system for the points GeoDataFrame
# GeoPandas doesn't do any transformations automatically when performing
# the spatial join. The layers are already in the same CRS (WGS84) so no
# transformation is needed.
points.crs = polys.crs
###print points.crs

#print nbrhoods
# Drop all columns except the name and polygon geometry

#boroughs = nbrhoods[["borough", "geometry"]]

# Perform the spatial join
result = geopandas.tools.sjoin(points, nbrhoods, how="left")
# Print the results...
``````

result: `ValueError: need at least one array to concatenate`. From researching this error, it seems that it happenned commonly when spatial join yields no results but it was supposed to be fixed by this version of GeoPandas. I assume it is fixed and it is another problem here.

### Attempt 3: GeoPandas 2 : (similar)

``````import pandas as pd
from shapely.geometry import Point
import geopandas as gp
from geopandas.tools import sjoin

#trip = pd.read_csv('TripRecordsReporttrips.csv', sep=',',header=None, usecols=[0, 4, 8, 9, 10, 11],names=['TripID', 'Date', 'StartLat', 'StartLon','EndLat','EndLon'])

geometry = [Point(-73.9868805930018, 40.7697167683218),
Point(-73.977673,40.722249)]
#geometry = [Point(xy) for xy in zip(trip['StartLon'], trip['StartLat'])]
#geometry2 = [Point(xy) for xy in zip(trip['EndLon'], trip['EndLat'])]
#trip = trip.drop(['StartLon', 'StartLat','EndLon','EndLat'], axis=1)
trip = ['broadway 150', 'crosby 10']
crs = {'init' :'epsg:4326'}
starts = gp.GeoDataFrame(trip, crs=FRC1.crs , geometry=geometry)
#ends  = gp.GeoDataFrame(trip, crs=crs, geometry=geometry2)

pointInPolys = sjoin(starts, FRC1, how='left',op="within")
``````

## Reproduction Materials

Shapefiles for neighborhoods is the same as here: https://geo.nyu.edu/catalog/nyu_2451_36743

Example input for points:

Start Lng................... Start Lat.................. End Lat ................... End Long................

-73.9446545764804 40.7796840845172 40.7693589967192 -73.9857926219702 -73.9574825763702 40.8116657733964 40.7834276940059 -73.9789965748787 -73.9574825763702 40.8116657733964 40.7839981278416 -73.9762181416154 -73.9938855171204 40.7355130703031 40.7112190324095 -73.9506489783525 -73.9634789898992 40.7588519932381 40.7544933950398 -73.9975795894861 -73.9907560497522 40.7364947075543 40.804403082633 -73.9372211694717 -73.9907560497522 40.7364947075543 40.804403082633 -73.9372211694717 -73.9939170330763 40.7464936766106 40.7191759823283 -73.9889881387353 -73.9907560497522 40.7364947075543 40.804403082633 -73.9372211694717 -73.9831633865833 40.718557218786 40.6432838990509 -73.7899511307478 -73.9887668564916 40.7186952022837 40.7759590522127 -73.9807621389627

*Points put into an array just for experiments - even through half of them are lat/long and half are long/lat, I'm trying both to see if the problem might be that the functions take it as one way and not the other.

[(-73.9446545764804, 40.7796840845172) (40.7693589967192 ,-73.9857926219702) (-73.9574825763702, 40.8116657733964 ) (40.7834276940059 ,-73.9789965748787 ) (-73.9574825763702, 40.8116657733964 ) (40.7839981278416 ,-73.9762181416154 ) (-73.9938855171204, 40.7355130703031 ) (40.7112190324095 ,-73.9506489783525 ) (-73.9634789898992, 40.7588519932381 ) (40.7544933950398 ,-73.9975795894861 ) (-73.9907560497522, 40.7364947075543 ) (40.804403082633 ,-73.9372211694717 ) (-73.9907560497522, 40.7364947075543 ) (40.804403082633 ,-73.9372211694717 ) (-73.9939170330763, 40.7464936766106 ) (40.7191759823283 ,-73.9889881387353 ) (-73.9907560497522, 40.7364947075543 ) (40.804403082633 ,-73.9372211694717 ) (-73.9831633865833, 40.718557218786 ) (40.6432838990509 ,-73.7899511307478 ) (-73.9887668564916, 40.7186952022837 ) (40.7759590522127, -73.9807621389627)]

It is simply a projection problem. The CRS of the taxi_zones.shp file is

``````c = fiona.open('taxi_zones.shp')
c.crs
{u'lon_0': -74, u'datum': u'NAD83', u'y_0': 0, u'no_defs': True, u'proj': u'lcc', u'x_0': 300000, u'units': u'us-ft', u'lat_2': 41.03333333333333, u'lat_1': 40.66666666666666, u'lat_0': 40.16666666666666}
``````

That corresponds to EPSG 2263 (NAD83 / New York Long Island (ftUS))

If I use the nyu_2451_36743.shp (nyu_2451_36743_WGS84.zip) shapefile (in WGS84 projection)

``````c = fiona.open('nyu_2451_36743.shp')
c.crs
{u'no_defs': True, u'datum': u'WGS84', u'proj': u'longlat'}
``````

Then

`````` NY_nbr_shpfile ='nyu_2451_36743.shp'
....
print coor_to_nbr(-73.9868805930018, 40.7697167683218, NY_nbr_shpfile)
[141]
OrderedDict([(u'OBJECTID', 142), (u'Shape_Leng', 0.03817589423), (u'Shape_Area', 7.565379e-05), (u'zone', u'Lincoln Square East'), (u'LocationID', 142), (u'borough', u'Manhattan')])
``````

works

• Yes! Thank you so much, Gene! true, the problem WAS indeed the PROJECTION type! (between NAD83 and WGS84 projections). It worked when I simply used the different versions of the files that are in the WGS84 projection, and also worked when I converted the projection to WGS84 within GeoPandas. Aug 1, 2017 at 21:23