Given list of points (lat, long), how to find all points within radius of a given point [closed]

Question

It's my first time working with GIS. I really wish I could share some code, but I don't even know how and where to begin. I have a Pandas dataframe with around 10,000 records. Each of them have a latitude and longitude. Assume that after applying some filters, I end up with 10 records relevant to me.

Now, I want to find all the records which have coordinates lying within 20 miles (radius...?) of ANY of these 10 records.

The dumbest brute force solution I thought of was using geopy and looping over the records and finding the distance. Then I came across this solution to use kdtree - https://stackoverflow.com/a/41155458. But the comment below it says that it isn't accurate.

And I also tried to look at geopandas, and use the buffer and unary_union functions - http://geopandas.org/geometric_manipulations.html. The only thing I have with me right now is the dataframe with the latitudes and longitudes.

How should I proceed?

Answer 1

1. Create GeoPandas geodataframes:

   import geopandas as gpd
   import shapely
   df # your pandas dataframe with 10k records
   df_filt # your filtered dataframe
   
   # Create geometries from your lat-lons
   geom_list = [shapely.geometry.Point(lon,lat) for lon,lat in zip(df["longitude" ,df["latitude"])] # check the ordering of lon/lat
   
   # create geopandas geodataframe
   gdf = gpd.GeoDataFrame(df, geometry=geom_list, crs={"init":"EPSG:4326"})
   

Do the same for your filtered data. If you're doing filtering in pandas, you can also do it on your geodataframe instead.

2. Reproject to an appropriate projected coordinate system (units in meters/miles, etc, not degrees!) depending on the location and extent of your records. I cannot answer that, but if you figure out an appropriate system and EPSG code, you can convert your gdf as follows:

   gdf_proj = gdf.to_crs({"init": "EPSG:<put code here>"})
   

3. Then buffer your geodataframe of the 10 records, use unary_union, then intersect your 10k record gdf with the union result. See how you go with those steps and put it as an answer if you figure it out.

Answer 2

Following the answer given by @sbphd, this is what I coded.

import geopandas as gpd

# with columns "id", "latitude", "longitude" - 10k records
df

gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(
        df["longitude"],
        df["latitude"],
    ),
    crs={"init":"EPSG:4326"},
)

# 10 records
filtered_df

filtered_gdf = gpd.GeoDataFrame(
    filtered_df, 
    geometry=gpd.points_from_xy(
        filtered_df["longitude"],
        filtered_df["latitude"],
    ),
    crs={"init":"EPSG:4326"},
)

# EPSG:3857 converts it to meters, correct?

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

# so 100 miles would be 160934 meters

x = filtered_gdf_proj.buffer(160934).unary_union

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

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