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

Question

我有一个df像这样的 Pandas 数据框：

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

我还有一个带有poly一些多边形的 Geopandas 数据框。现在，我只想绘制某个多边形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()

给了我这个图像。

请注意，您可以以其他更快的方式制作布尔掩码，例如纬度是否高于多边形的最高点。这些本身可能并不完美，但可以减少问题，这样您就不必intersects()多次拨打电话。

[编辑：如果您的多边形是矩形，]另一种方法（正如您在问题中所建议的那样）是仅“裁剪”边界多边形周围的图像。这是一个更快的解决方案，因为您不必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])) )

给出以下内容：