CODEWITHSUNDEEP
pythonpandasgeometryshapelygeopandas
Tashay Green
Tashay Green

Reputation: 105

How to Find a Point within a Polygon?

I am trying to find a point within polygons of a shapefile.

I need to write a loop that can loop over the polygons and return the index of the polygon in which the point is located.

How would I write a loop to find out which polygon the point is in?

Here's what I have written so far: 

import pandas as pd
import pylab as pl 
import os
import zipfile
import geopandas as gp
import shapely

%pylab inline

# Read in the shapefile 
ct_shape = gp.read_file(path)

# Segmented the file so it only contains Brooklyn data & set projection
ct_latlon = ct_shape[ct_shape.BoroName == 'Brooklyn']
ct_latlon = ct_latlon.to_crs({'init': 'epsg:4326'}) 
ct_latlon.head()

# Dataframe image
[Head of the dataframe image][1]: https://i.sstatic.net/xAl6m.png


# Created a point that I need to look for within the shapefile
CUSP = shapely.geometry.Point(40.693217, -73.986403)

The output could be something like this: '3001100' (the BCTCB2010 of the correct polygon)

Upvotes: 1

Views: 4297

Answers (3)

user7287085
user7287085

Reputation:

Something that might be of interest, in addition to your accepted answer: you can also take advantage of geopandas's built-in Rtree spatial indexing for fast intersection/within queries.

spatial_index = gdf.sindex
possible_matches_index = list(spatial_index.intersection(polygon.bounds))
possible_matches = gdf.iloc[possible_matches_index]
precise_matches = possible_matches[possible_matches.intersects(polygon)]

From this tutorial. The example returns which points intersect a single polygon, but you can easily adapt it to your example of a single point with multiple polygons.

Upvotes: 0

Tashay Green
Tashay Green

Reputation: 105

I solved it in one line of code. No loop necessary.

Posting for anyone else that may be interested: 

# Setting the coordinates for the point
CUSP = shapely.geometry.Point((-73.986403, 40.693217,)) # Longitude & Latitude

 # Printing a list of the coords to ensure iterable 
 list(CUSP.coords)

 # Searching for the geometry that intersects the point. Returning the index for the appropriate polygon. 
 index = ct_latlon[ct_latlon.geometry.intersects(CUSP)].BCTCB2010.values[0]

Upvotes: 1

arthur
arthur

Reputation: 2399

I would use a GeoDataFrame sjoin.

Below is a short example, I have one city corresponding to Paris coordinates and I will try to match it with a country that is in the countries GeoDataFrame.

import geopandas as gpd
from shapely.geometry.point import Point

# load a countries GeoDataFrame given in GeoPandas
countries = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\
                .rename(columns={"name":"country_name"})

#making a GeoDataFrame with your city
paris = Point( 2.35, 48.85)
cities = gpd.GeoDataFrame([{"city" : "Paris", "geometry":paris} ])

In [33]: cities  
Out[33]: 
    city            geometry
0  Paris  POINT (2.35 48.85)

#now we left_join cities and countries GeoDataFrames with the operator "within" 

merging = gpd.sjoin(cities, countries, how="left", op="within")

In [34]: merging 
Out[34]: 
    city            geometry  index_right continent  gdp_md_est iso_a3  \
0  Paris  POINT (2.35 48.85)           55    Europe   2128000.0    FRA   

  country_name     pop_est  
0       France  64057792.0

We see that the paris Point has been found inside the polygon of the country at the index 55 in the countries GeoDataFrame which is France : 

In [32]: countries.loc[55]
Out[32]: 
continent                                                  Europe
gdp_md_est                                              2.128e+06
geometry        (POLYGON ((-52.55642473001839 2.50470530843705...
iso_a3                                                        FRA
country_name                                               France
pop_est                                               6.40578e+07
Name: 55, dtype: object

Thus, if you have a list of points instead of just one, you just have to create a bigger cities GeoDataFrame.

Upvotes: 0

Related Questions