4

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

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

poly我还有一个带有一些多边形的 Geopandas 数据框。现在，我只想绘制某个多边形df内的点。所以我应该能够做类似的事情，Shapelypoly.intersects(p)在哪里。但我做错了什么；pPoint

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

传递带有纬度/经度点的数据框并将它们叠加到的最佳方法是poly什么？请注意，我可以定义一个最小/最大纬度/经度范围，但这也会poly在（更大）边界框的外部但内部打印点。

4

2 回答 2

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

给出以下内容：

在此处输入图像描述

于 2016-09-15T01:16:55.870 回答
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:30:36.913 回答