如何删除 X 和 Y 坐标在多边形之外的数据框行
How to drop dataframe rows where X and Y coordinates are outside of polygon
我正在尝试解决以下问题。假设一个数据框（从 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)
这很好用。但是我不知道如何继续这个。
非常感谢帮助
----编辑----
请参阅下面的更新解决方案
我对shapely不是很熟悉。也许他们有真正的 pandas 支持。 Afaik，他们支持向量化的 numpy 函数，所以我不会感到惊讶。
找出哪些点在给定多边形内的一种方法是使用 pandas apply() 函数：
import pandas as pd
from shapely.geometry import Polygon, Point
#your dataframe of points
df = pd.DataFrame([[0, 0, 0], [1, 2, 3], [2, 2, 2], [3, 2, 1] ], columns = list("XYZ"))
#your polygon points
polygon1_list = [(1, 1), (1, 3), (3, 3), (3, 1)]
#adding a column that contains a boolean variable for each point
df["polygon1"] = df.apply(lambda row: Polygon(polygon1_list).contains(Point(row["X"], row["Y"])), axis = 1)
print(df)
我的玩具数据集的输出
X Y Z polygon1
0 0 0 0 False
1 1 2 3 False
2 2 2 2 True
3 3 2 1 False
在 shapely 中，contains 真正意味着在多边形内，这不包括边界。如果你想包括边框，你应该使用 intersects
df["polygon1"] = df.apply(lambda row: Polygon(polygon1_list).intersects(Point(row["X"], row["Y"])), axis = 1)
现在你的问题的答案很简单了。只需将包含 False 的行放到这个新列中：
df = df.drop(df[~df["polygon1"]].index)
不幸的是，您仍然需要遍历多边形列表。如果有人知道如何在没有（显式）循环的情况下测试所有点和所有多边形，那将会很有趣。我见过一个 MultiPolygon 构造函数 class on their website，所以也许将所有多边形组合成一个 class 就可以了。但请提前测试这是一个有效的选择。如果 MultiPolygon 的成员沿着一条直线接触无限多的点，则该 MultiPolygon 无效。
编辑：看起来，在 Python 2.7 中这不起作用。
@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函数快多了。虽然公平地说 non-geo pandas 解决方案还必须将 X 和 Y 转换为形状点元素，然后执行交集评估
我无法模仿 Python 2.7 中建议的 。所以这是我必须做出的细微差别才能让它在 Python 2.7.
中工作
from shaply.geometry.polygon import Polygon
inside = Polygon(poly_points).contains_points(zip(df.X.values, df.Y.values))
df['inside'] = inside
df = df.drop(df[~df['inside']].index)
看来老版本contains_points单点运行有问题。所以我将其设置为读取所有要点并将该列表附加为新列。
我正在尝试解决以下问题。假设一个数据框（从 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)
这很好用。但是我不知道如何继续这个。
非常感谢帮助
----编辑---- 请参阅下面的更新解决方案
我对shapely不是很熟悉。也许他们有真正的 pandas 支持。 Afaik，他们支持向量化的 numpy 函数，所以我不会感到惊讶。
找出哪些点在给定多边形内的一种方法是使用 pandas apply() 函数：
import pandas as pd
from shapely.geometry import Polygon, Point
#your dataframe of points
df = pd.DataFrame([[0, 0, 0], [1, 2, 3], [2, 2, 2], [3, 2, 1] ], columns = list("XYZ"))
#your polygon points
polygon1_list = [(1, 1), (1, 3), (3, 3), (3, 1)]
#adding a column that contains a boolean variable for each point
df["polygon1"] = df.apply(lambda row: Polygon(polygon1_list).contains(Point(row["X"], row["Y"])), axis = 1)
print(df)
我的玩具数据集的输出
X Y Z polygon1
0 0 0 0 False
1 1 2 3 False
2 2 2 2 True
3 3 2 1 False
在 shapely 中，contains 真正意味着在多边形内，这不包括边界。如果你想包括边框，你应该使用 intersects
df["polygon1"] = df.apply(lambda row: Polygon(polygon1_list).intersects(Point(row["X"], row["Y"])), axis = 1)
现在你的问题的答案很简单了。只需将包含 False 的行放到这个新列中：
df = df.drop(df[~df["polygon1"]].index)
不幸的是，您仍然需要遍历多边形列表。如果有人知道如何在没有（显式）循环的情况下测试所有点和所有多边形，那将会很有趣。我见过一个 MultiPolygon 构造函数 class on their website，所以也许将所有多边形组合成一个 class 就可以了。但请提前测试这是一个有效的选择。如果 MultiPolygon 的成员沿着一条直线接触无限多的点，则该 MultiPolygon 无效。
编辑：看起来，在 Python 2.7 中这不起作用。
@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函数快多了。虽然公平地说 non-geo pandas 解决方案还必须将 X 和 Y 转换为形状点元素，然后执行交集评估
我无法模仿 Python 2.7 中建议的 Python 2.7.
from shaply.geometry.polygon import Polygon
inside = Polygon(poly_points).contains_points(zip(df.X.values, df.Y.values))
df['inside'] = inside
df = df.drop(df[~df['inside']].index)
看来老版本contains_points单点运行有问题。所以我将其设置为读取所有要点并将该列表附加为新列。