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

from How to map a point to a polygon label from a shapefile:

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
        # 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

# Load the 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

FRC1  = gp.read_file(ny_nbr_shpfile)

#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),
#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)

#starts.head() ###print 

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')
{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')
{u'no_defs': True, u'datum': u'WGS84', u'proj': u'longlat'}


 NY_nbr_shpfile ='nyu_2451_36743.shp'
 print coor_to_nbr(-73.9868805930018, 40.7697167683218, NY_nbr_shpfile)
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')])


  • 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. – pentatomic Aug 1 '17 at 21:23

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.