2

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: enter image description here

I now want to use Python to essentially create a bit like the Select by Location tool in ArcGIS Pro https://pro.arcgis.com/en/pro-app/2.8/tool-reference/data-management/select-layer-by-location.htm

So for example the script would work like:

  • 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 this so far

import geopandas as gpd
import pandas as pd
import shapely

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).unary_union

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

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

But I need a way of defining which ID I want to set the 10 metre buffer from

3
  • 1
    Welcome to GIS SE. As a new user, please take the Tour. Programming questions here are expected to contain code. Please choose a code platform and make a coding attempt, then, if there's an issue, you can Edit this Question to focus on that issue.
    – Vince
    Commented Jun 16, 2022 at 12:48
  • @Vince see update Commented Jun 16, 2022 at 13:15
  • geopandas has a distance method to calculate the distance and you can sort it, apply threshold from there, and it returns a boolean array. have you used the distance method?
    – sutan
    Commented Jun 16, 2022 at 14:52

1 Answer 1

1

You can create a dataframe of one point, buffer the point and spatial join it to all other points. Spatial join is very fast.

import geopandas as gpd

df = gpd.read_file(r'/home/bera/GIS/Data/LMV-data_2021-10/ok_riks_Sweref_99_TM_shape/oversikt/riks/gs_riks.shp')

onepoint = df.loc[df['id']==10].copy() #Select one point
onepoint.geometry = onepoint.geometry.buffer(10000)

result = sorted(onepoint.sjoin(df.loc[df['id']!=10]).id_right.tolist())

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.