I have created a test shapefile containing 15 point features in EPSG:2157 and exported it geojson. Each point has been assigned an ID - e.g. 1, 2 ,3 , etc. They look like so:
{
"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
I now want to use Python to essentially:
- Specify the ID of the point of interest
- Add a search distance in metres
- Print the ID's of the points within the specified distance and their total distance from the point of interest
I have tried geopandas
so far to get me going as per https://gis.stackexchange.com/questions/349637/given-list-of-points-lat-long-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])
But this is just printing my original geopandas dataframe with all 15 IDs and longitude/latitudes,
I need a way of specifying which ID I want from the dataframe, set the 10 metre buffer on it and from that print whichever of the remaining 14 points ID and distance from that point.
How do I go about such a thing?
CodePudding user response:
This approach first calculates the distance from the selected point and then filters to the search distance:
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)
Output:
id geometry distance
0 1 POINT (-657346.500 7286537.676) 0.000000
1 2 POINT (-657339.334 7286538.871) 7.264817