在一个数据框中围绕不同数据框（例如酒店）中的点查找实体（例如餐馆）的数量（坐标计数问题）

Question

对于一个项目，我们正在尝试计算（并命名）一个数据框中的点数，这些点在另一个数据框中的点周围，具有给定的特定半径。我们尝试了很多，但通过手动计算 Tableau 中的点来验证我们的解决方案并没有达到令人满意的解决方案。虽然我们相当接近。 我们有两个数据框。 一个数据框有大约 70k 行和 50 列，具有唯一的酒店 ID、纬度、经度、名称和酒店的不同信息（例如“has_desk”是/否等）。 另一个有大约 25,000 行和 9 列，具有唯一的机构 ID、纬度、经度、名称、设施类型（例如“餐厅”与“酒吧”）和其他信息，例如美食和 vegan_available 等

由于数据集的大小，为每家酒店计算到每家餐厅等的距离的嵌套循环似乎是不可能的。 出于计算原因，使用六边形而不是酒店周围的真实圆圈似乎也是个好主意。

输入：

radius_in_m = 200

df_hotels:

    id  lat lon name
0   1   50.600840   -1.194608   Downtown Hotel
1   2   50.602031   -10.193503  Hotel 2
2   3   50.599579   -10.196028  Hotel 3

df_poi:

    id  lat         lon         name                    amenity
0   451152  51.600840   -0.194608   King of Prussia         restaurant
1   451153  51.602031   -0.193503   Central Restaurant  restaurant
2   451154  51.599579   -0.196028   The Catcher in the Rye  bar

想要的结果： df_hotels_new:

    id  lat lon name      num_restaurants       restaurants_list     num_bar     bars_list  
0   1   50.600840   -1.194608   Downtown Hotel        2         [451152, 451153]      0       []
1   2   50.602031   -10.193503  Hotel 2               0         []                    1       [451154]
2   3   50.599579   -10.196028  Hotel 3               0         []                    0       []

在示例中，前两家餐厅距离第一家酒店不到 200 米。该计数被添加到新列中。包含计算在内的两家餐厅 ID 的列表将添加到另一列。第三个是酒吧，因此不计入餐厅。请注意，示例中的纬度/经度完全是虚构的，实际上不在 200 米半径内。

迄今为止最成功的尝试是这个，但它大多高估了餐厅的数量。它也没有在另一列中列出餐厅/酒吧/等，但我们已经让它起作用了。通过它，我们能够看到半径似乎比指定的“稍微”（大约 1.5 倍）大，并且可能也有一点偏移。这可能是四舍五入或地图投影错误吗？

import geopandas as gpd
from shapely.geometry import Point
from shapely.ops import transform
from functools import partial
import pyproj
import math

# Define the conversion factor from meters to degrees based on the latitude
def meters_to_degrees(meters, latitude):
    proj_meters = pyproj.CRS("EPSG:3857")  # meters
    proj_latlon = pyproj.CRS("EPSG:4326")  # degrees
    transformer = pyproj.Transformer.from_crs(proj_meters, proj_latlon, always_xy=True)
    lon, lat = transformer.transform(meters, 0)
    lat_dist_per_deg = 111132.954 - 559.822 * math.cos(2 * math.radians(latitude)) + 1.175 * math.cos(4 * math.radians(latitude))
    lon_dist_per_deg = 111412.84 * math.cos(math.radians(latitude))
    lat_degrees = meters / lat_dist_per_deg
    lon_degrees = meters / lon_dist_per_deg
    return lat_degrees, lon_degrees




# Convert the hotels DataFrame to a GeoDataFrame with a Point geometry column
hotels_geo = gpd.GeoDataFrame(df_hotels, geometry=gpd.points_from_xy(df_hotels["longitude"], df_hotels["latitude"]))

# Convert the poi/restaurant DataFrame to a GeoDataFrame with a Point geometry column
poi_geo = gpd.GeoDataFrame(df_poi, geometry=gpd.points_from_xy(df_poi["longitude"], df_poi["latitude"]))

# Create an R-tree spatial index for the df_poi GeoDataFrame
df_poi_sindex = poi_geo.sindex

# Define the radius of the search in meters
radius_meters = 200

# Loop through each row in hotels_geo
for index, row in hotels_geo.iterrows():
    # Convert the radius from meters to degrees based on the latitude
    lat, lon = row["latitude"], row["longitude"]
    lat_deg, lon_deg = meters_to_degrees(radius_meters, lat)
    
    # Use the R-tree spatial index to find the df_poi rows within the search radius
    candidate_indices = list(df_poi_sindex.intersection(row.geometry.buffer(lon_deg).bounds))

    # Filter the street_test rows to only those within the search radius
    candidate_rows = poi_geo.iloc[candidate_indices]

    # Group the candidate rows by amenity and count the occurrences
    counts = candidate_rows.groupby("amenity").size().to_dict()

    # Add the counts as new columns in the df_hotels DataFrame
    for amenity_type, count in counts.items():
        df_hotels.at[index, amenity_type] = count

    # Print progress
    if index % 10000 == 0:
        print(f"Processed {index} rows")

# Replace NaN values with 0
airbnb_test.fillna(value=0, inplace=True)