地理数据框以选择某些ID并在指定距离内找到几何形状

Question

我创建了一个测试形状文件，其中包含EPSG：2157中的15个点功能，并导出了Geojson。每个点都被分配了一个ID-例如1、2、3等。它们看起来像：

{
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {
                "id": "1"
            },
            "geometry": {
                "type": "Point",
                "coordinates": [
                    -5.905044078826904,
                    54.609987802465916
                ]
            }
        },
      {
            "type": "Feature",
            "properties": {
                "id": "11"
            },
            "geometry": {
                "type": "Point",
                "coordinates": [
                    -5.902683734893799,
                    54.60972062159888
                ]
            }
        }
    ]
}

etc

我现在想使用python来本质上：

• 指定感兴趣点的ID
• 添加搜索距离，以仪表
• 打印iD指定距离内的点及其距离感兴趣点的总距离，

我尝试了geopandas到目前为止，以使我按照 https://gis.stac.stackexchange.com/questions/questions/349637/given - 列表lat-lat-long-how-how-to-find-all-points-within-radius of-a-give ，

import geopandas as gpd
import pandas as pd

input_file = 'C:/test/points.geojson'
df = gpd.read_file(input_file)
df['lon'] = df['geometry'].x
df['lat'] = df['geometry'].y


gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(
        df["lon"],
        df["lat"],
    ),
    crs={"init":"EPSG:2157"},
)

print(gdf)

gdf_proj = gdf.to_crs({"init": "EPSG:3857"})


x = gdf_proj.buffer(10)

neighbours = gdf_proj["geometry"].intersection(x)

# print all the nearby points
print(gdf_proj[~neighbours.is_empty])

但这只是在打印我的原始geopandas dataframe，其中所有15个ID和经度/经度/纬度，

我需要一种从数据框架中指定所需的ID的方法，在其上设置了10米缓冲区，以及从该点剩下的14点ID和该点的距离中的任何一个打印中。

我该如何做这样的事情？

Answer 1

该方法首先计算从选定点，然后过滤到搜索距离的距离：

import geopandas as gpd

input_file = 'test.geojson'
gdf = gpd.read_file(input_file).to_crs('EPSG:3857')
my_id = 1
search_distance = 10
selected_feature = gdf.loc[gdf['id'].astype(int) == my_id]
result = (
    gdf
    .assign(distance=gdf.apply(lambda x: x.geometry.distance(selected_feature.geometry.iloc[0]), axis=1))
    .query(f'distance <= {search_distance}')
)
print(result)

输出：

  id                         geometry  distance
0  1  POINT (-657346.500 7286537.676)  0.000000
1  2  POINT (-657339.334 7286538.871)  7.264817