如何将pandas中的lat/lon个点进行转换,看它们是否落在某些边界多边形中?

How to convert lat/lon points in pandas and see whether they fall in some boundary polygons?

我有一个 Pandas 数据框 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))

传递具有 lat/lon 点的数据帧并将它们叠加到 poly 上的最佳方法是什么?请注意,我可以定义 min/max lat/lon 的范围,但这也会在 poly 之外但在(更大的)边界框内打印点。

你的起点:

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

[编辑:如果您的多边形是矩形,] 另一种方式(如您在问题中所建议的那样)就是 'crop' 边界多边形周围的图像。这是一个更快的解决方案,因为您不必一遍又一遍地调用 intersects() 函数。要 trim 基于边界多边形的图像,您可以在 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])) )

给出以下内容:

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

然后它以不同的颜色绘制多边形及其内部和外部的点,就像您想要的那样。