如何转换大熊猫经纬度点并查看它们是否落在某些边界多边形中？

Question

我有一个熊猫数据帧df这样的：

id lat lon 
jhg 2.7 3.5 
ytr 3.1 3.5 
... 

我也有一个Geopandas数据帧poly一些多边形。现在，我想只绘制点df中的内部的一些多边形。所以我应该可以做点像poly.intersects(p)，其中p是Shapely Point。但我做错了什么;

from shapely.geometry import Point 
for index, row in df.iterrows(): 
    t = poly.intersects(Point(row.lon, row.lat)) 

什么是通过与经/纬度点的数据框并绘制它们叠加到poly最好的方法是什么？请注意，我可以定义一个最小/最大纬度/经度范围，但也可以打印poly以外的点，但位于（较大）边界框内。

Answer 1

你的出发点：

import pandas as pd 
from shapely.geometry import box 
import matplotlib.pyplot as plt 

from matplotlib.collections import PatchCollection 
from matplotlib.patches import Polygon 
from shapely.geometry import Point 
import seaborn as sns 
import numpy as np 

# some pretend data 
data = {'lat':[2.7,3.5,1.4,2.3,.9,1.9], 'lon':[1.2,.9,1.9,2.2,3,1.1]} 
df = pd.DataFrame(data) 

# the 'bounding' polygon 
poly = box(1,1,2,2) 
patches = PatchCollection([Polygon(poly.exterior)], facecolor='red', linewidth=.5, alpha=.5) 


# plot the bounding box 
fig, ax = sns.plt.subplots(1, figsize=(4,4)) 
ax.add_collection(patches, autolim=True) 

# plot the lat/lon points 
df.plot(x='lat',y='lon', kind='scatter',ax=ax) 
plt.show() 

的数字看起来是这样的：

一种方式来摆脱不必要的点是使用一个布尔面膜：

#probably more efficient ways to do this, but this works 
mask = [poly.intersects(Point(lat,lon)) for lat,lon in zip(df.lat,df.lon)] 
df = df[mask] 

# make new plot (must make a new 'patch' object) 
patches1 = PatchCollection([Polygon(poly.exterior)], facecolor='red', linewidth=.5, alpha=.5) 
fig1, ax1 = sns.plt.subplots(1, figsize=(4,4)) 
ax1.add_collection(patches1, autolim=True) 

# make the axis bounds the same 
ax1.set_xlim(ax.get_xlim()) 
ax1.set_ylim(ax.get_ylim()) 

# plot the lat/lon points 
df.plot(x='lat',y='lon', kind='scatter',ax=ax1) 
plt.show() 

将此图像给我。

请注意，你可以做其他，更快的方式一个布尔面具等是否LAT是上述多边形中的最高点。这些可能并不完美，但可以减少问题，因此您没有多次致电intersects()。

[编辑：如果你的多边形是一个矩形，]另一种方式（就像你在你的问题中所建议的那样）将只是围绕边界多边形“裁剪”图像。这是一个更快的解决方案，因为您不必一次又一次调用该函数。要修剪基于边界多边形的图像，你可以plt.plot()之前插入这一权利：

ax.set_xlim((np.min(poly.exterior.xy[0]),np.max(poly.exterior.xy[0]))) 
ax.set_ylim((np.min(poly.exterior.xy[1]),np.max(poly.exterior.xy[1]))) 

提供了以下：

Answer 2

This tutorial似乎做你想要什么。它还利用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)] 

然后它绘制多边形和它内部的点和外部它以不同的颜色，就像你想要做的一样。