0

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?

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.

  1. 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>"})
    
  2. 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.

  • Added an answer with some code. Could please you take a look? Thanks! – Jeet Parekh Feb 5 at 19:31
0

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])
  • Looks good, but does it work? I didn't know about '''gpd.points_from_xy''' - that's a much nicer way to create the geometries. I'm not sure of the accuracy of EPSG:3857 in this context. I know enough to know that I don't fully understand projections... but assuming your 10 mile radius doesn't need to be perfectly accurate, you'll probably be fine. Otherwise, depending on the area you cover, use the local UTM zone or country grid. Usually you reference data frames with '''.loc''' or '''.iloc''', but your solution probably works too. – sbphd Feb 6 at 8:57

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