1

我正在尝试在 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）

4

3 回答 3

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]
于 2016-11-28T17:42:39.063 回答
0

我会使用 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已在国家的多边形内找到，位于countriesGeoDataFrame 中的索引 55 处，即 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

因此，如果你有一个点列表而不是一个点，你只需要创建一个更大的citiesGeoDataFrame。

于 2016-11-28T17:24:26.543 回答
0

除了您接受的答案之外，可能会感兴趣的东西：您还可以利用 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)]

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

于 2016-12-12T20:26:14.400 回答