python-3.x - 如何删除 X 和 Y 坐标位于多边形之外的数据框行

Question

我正在尝试解决以下问题。让我们假设一个数据帧(从 txt 文件加载)具有以下结构(和数千行):

foo.head()

         X            Y       Z 
 0  125417.5112  536361.8752 -1750.0
 1  127517.7647  533925.8644 -1750.0
 2  128144.1000  533199.4000 -1750.0
 3  128578.8385  532904.9288 -1750.0
 4  125417.5112  536361.8752 -1750.0
 ....

数据代表 X、Y 和 Z 坐标。

我还有一组定义闭合多边形的点。这些位于 numpy 数组中:

polypoints

array([[ 125417.5112,  536361.8752],
       [ 127517.7647,  533925.8644],
       [ 128144.1   ,  533199.4   ],
       ....
       [ 125417.5112,  536361.8752]])

如何过滤数据框以删除不属于闭合多边形内的行？

我尝试使用shapely.geometry polygon定义多边形。通过这样做:

poly = Polygon(polypoints)

这很好用。但我不知道如何继续下去。

非常感谢您的帮助

----编辑----请参阅下面的更新解决方案

Answer 1

@MrT 建议的原始解决方案效果很好。然而，根据 @Rutger Kassies 的建议查看 geopandas，我还找到了另一个解决方案。第一个需要安装 geopandas 包。那么下面的代码对我有用:

import geopandas as gpd
from shapely.geometry import Point, Polygon, MultiPolygon
# load the data that should be cropped by the polygon
# this assumes that the csv file already includes 
# a geometry column with point data as performed below
dat_gpd = gpd.GeoDataFrame.from_csv(r'{}\data_to_crop.csv'.format(savedir), sep='\t')

# load the data of the polygon as a dataframe
arr_df = pd.DataFrame(data, columns=['X','Y','Z'])

# make shapely points out of the X and Y coordinates
point_data = [Point(xy) for xy in zip(arr_df.X, arr_df.Y)]

# assign shapely points as geometry to a geodataframe
# Like this you can also inspect the individual points if needed
arr_gpd = gpd.GeoDataFrame(arr_df, geometry=point_data)

# define a shapely polygon from X and Y coordinates of the shapely points
polygo = Polygon([[p.x, p.y] for p in arr_gpd.geometry])

# assing defined polygon to a new dataframe
pol_gpd= gpd.GeoDataFrame()
pol_gpd['geometry'] = None
pol_gpd.loc[0,'geometry'] = polygo

# define a new dataframe from the spatial join of the dataframe with the data to be cropped
# and the dataframe with the polygon data, using the within function.
dat_fin = gpd.sjoin(dat_gpd, pol_gpd, op = 'within')

希望这对遇到类似问题的人有所帮助。此外，可以找到有关空间连接的更多信息 on the geopandas website 。请注意，此功能不需要多边形之间的操作，但也适用于点和多边形

--编辑--

%timeit gpd.sjoin(dat_gpd, pol_gpd, op = 'within')
31.8 s ± 108 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit dat_gpd['inpoly'] = dat_gpd.apply(lambda row: polygo.intersects(Point(row["X"], row["Y"])), axis = 1)
1min 26s ± 389 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

看来geo-pandas功能要快得多。虽然公平地说，非地理 Pandas 解决方案还必须将 X 和 Y 转换为形状点元素，然后执行相交评估