如何在多边形内找到一个点？答案

【问题标题】：How to Find a Point within a Polygon?如何在多边形内找到一个点？
【发布时间】：2018-12-12 20:28:51
【问题描述】：

我正在尝试在 shapefile 的多边形内找到一个点。

我需要编写一个循环来遍历多边形并返回该点所在多边形的索引。

我将如何编写一个循环来找出该点所在的多边形？

这是我目前所写的：

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.stack.imgur.com/xAl6m.png


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

输出可能是这样的：'3001100'（正确多边形的 BCTCB2010）

【问题讨论】：

有什么问题？我在您的帖子中没有看到任何问题，也没有看到您遇到的问题的描述。
@OliverW。我需要编写一个循环来找到多边形内点的索引。
假设总共有 N 条边。您是否可以使用轮流尝试每个多边形的解决方案来解决总工作量 O(N) 还是您想要更高效，例如 O(log(N) + K)，其中 K 是穿过水平线的边数? （经过一些预处理。）

标签： python pandas geometry shapely geopandas

【解决方案1】：

我用一行代码解决了它。不需要循环。

为其他可能感兴趣的人发帖：

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

【讨论】：

【解决方案2】：

我会使用 GeoDataFrame sjoin。

下面是一个简短的示例，我有一个与巴黎坐标相对应的城市，我将尝试将其与国家 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

我们看到巴黎Point 已在法国countries GeoDataFrame 中索引55 的国家多边形内找到：

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

因此，如果您有一个点列表而不仅仅是一个点，您只需创建一个更大的cities GeoDataFrame。

【讨论】：

【解决方案3】：

除了您接受的答案之外，您可能会感兴趣：您还可以利用 geopandas 的内置 Rtree 空间索引来进行快速交叉/组内查询。

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)]

来自this tutorial。该示例返回哪些点与单个多边形相交，但您可以轻松地将其调整为具有多个多边形的单个点的示例。

【讨论】：