400万行数据：如何快速计算每行30km半径内的样本点数？

针对你要计算每个样本30km半径内数量的需求，以下是几个实用的优化方向，比单纯调整网格更高效：

一、网格策略的优化（而非盲目细化） #

你原有的奇偶网格思路是为了避免边界遗漏，但可以调整得更合理：

• 把网格单元设为略小于30km的尺寸（比如25km），换算成经纬度的话：法国纬度范围42°-51°，1°纬度≈111km，所以纬度间隔取0.22°（≈25km）；1°经度≈111×cos(纬度)km，法国区域经度间隔取0.35°-0.4°（≈25km）。
• 不用搞两套奇偶网格，直接给每个样本标记所属网格的坐标（比如(x_grid, y_grid)），然后对每个样本，只需要计算自身网格+相邻8个网格内的样本距离——因为30km范围最多覆盖9个25km的网格，这样能直接排除90%以上的无关样本，大幅减少计算量。
• 别把网格弄太小，比如小于10km，会导致网格数量暴增，反而增加管理和查询的开销。

二、用空间索引替代网格遍历（最推荐） #

空间索引是处理大规模地理数据的核心优化手段，比网格方法效率高得多：

• 先转平面坐标：经纬度是球面坐标，直接算距离误差大且慢，把所有样本转成法国对应的UTM投影（比如UTM 31N或32N，覆盖法国大部分区域），转成米为单位的平面坐标，这样可以用欧氏距离计算，速度更快。
• 构建KDTree/R-tree索引：以Python为例，用scipy.spatial.KDTree或者geopandas的空间索引：

  import numpy as np
  from scipy.spatial import KDTree
  import geopandas as gpd
  
  # 假设df是你的数据，包含lon, lat列
  gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat))
  # 转UTM投影（法国用EPSG:32631或32632）
  gdf_utm = gdf.to_crs(epsg=32631)
  # 提取平面坐标数组
  coords = np.array(list(zip(gdf_utm.geometry.x, gdf_utm.geometry.y)))
  # 构建KDTree
  tree = KDTree(coords)
  # 查询每个点30km（30000米）内的点数，减去自身（因为包含自己）
  gdf['density_count'] = tree.query_ball_point(coords, r=30000, return_length=True) - 1
  

  这种方法的时间复杂度是O(n log n)，比两两计算的O(n²)快几个数量级，400万条数据完全能处理。
• 如果用数据库的话，直接用PostGIS的ST_DWithin函数，配合空间索引，能高效完成批量查询。

三、计算前的粗过滤 #

如果不想转投影，用球面距离计算的话，可以先做粗过滤：

• 对每个样本，先计算经纬度的最大允许差值：纬度差≤30/111≈0.27°，经度差≤30/(111×cos(lat))≈0.43°（法国纬度对应的cos值约0.62-0.75）。
• 先过滤掉经纬度差超过这个范围的样本，再用Haversine公式精确计算距离，这样能减少需要精确计算的样本对数量。

四、其他小技巧 #

• 矢量化计算：避免用Python循环逐行处理，尽量用numpy、pandas的矢量化操作，利用底层C实现加速。
• 分块处理：如果内存不够，把数据按网格或区域分块，每次处理一块内的样本和相邻块的样本，分批次计算后合并结果。

内容的提问来源于stack exchange，提问作者Anthony Calonne