如何在多边形中找到一个点?

How to Find a Point within a Polygon?

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

我需要编写一个循环来遍历多边形和 return 点所在多边形的索引。

我将如何编写一个循环来找出该点位于哪个多边形?

以下是我到目前为止所写的内容:

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)

我会使用 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。

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

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

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

除了您接受的答案之外,您可能还感兴趣:您还可以利用 geopandas 的内置 Rtree 空间索引进行快速 intersection/within 查询。

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。该示例 returns 点与单个多边形相交,但您可以轻松地将其调整为具有多个多边形的单个点示例。