使用大规模数据集的测地线和极坐标有效计算两点之间的距离

Question

我创建了一个大约 10_000 行的虚拟 sample，其中有两列（pikcup 和下车位置点 - 编码为字符串）。

我使用以下命令将样本读入极坐标数据帧：

df = pl.read_csv("./taxi_coordinates.csv")

我想使用模块有效地计算这些点之间的距离

from geopy.distance import geodesic

请注意，我正在尝试寻找最有效的方法，因为我的原始样本超过 3000 万行。

我的方法使用

map_rows()

def compute_coordinates_v2(df:pl.DataFrame, col:str) -> pl.DataFrame:
    target_col:str = 'pu_polygon_centroid' if col == 'pickup' else 'do_polygon_centroid'
    location_data:str = f'{col}_location_cleaned'
    coordinates:str = f'{col}_coordinates'
    df = df.with_columns(
        pl.col(target_col).str.replace_all(r'POINT \(|\)', '').alias(location_data)
    ).with_columns(
        pl.col(location_data).str.split(' ').alias(coordinates)
    )
    return df

df = compute_coordinates_v2(df, 'pickup')
df = compute_coordinates_v2(df, 'dropoff')

以上操作会生成两列列表类型

shape: (5, 1)
┌───────────────────────────────────┐
│ pickup_coordinates                │
│ ---                               │
│ list[str]                         │
╞═══════════════════════════════════╡
│ ["-73.95701169835736", "40.78043… │
│ ["-73.95701169835736", "40.78043… │
│ ["-73.95701169835736", "40.78043… │
│ ["-73.9656345353807", "40.768615… │
│ ["-73.9924375369761", "40.748497… │
└───────────────────────────────────┘
shape: (5, 1)
┌───────────────────────────────────┐
│ dropoff_coordinates               │
│ ---                               │
│ list[str]                         │
╞═══════════════════════════════════╡
│ ["-73.9656345353807", "40.768615… │
│ ["-73.95701169835736", "40.78043… │
│ ["-73.95701169835736", "40.78043… │
│ ["-73.9924375369761", "40.748497… │
│ ["-74.007879708664", "40.7177727… │
└───────────────────────────────────┘

现在要计算距离，我使用以下代码

func

def compute_centroid_distance_v2(row):
    if (row[0][0]) and (row[0][1]) and (row[1][0]) and (row[1][1]):
        centroid_distance = geodesic(
            (row[0][1], row[0][0]), #(latitude, longitude)
            (row[1][1], row[1][0])
        ).kilometers
    else:
        centroid_distance = 0.0
    return centroid_distance

df = df.with_columns(
        df.select(["pickup_coordinates", "dropoff_coordinates"]).map_rows(compute_centroid_distance_v2).rename({'map': "centroid_distance"})
    )

以 3000 万行为基准，

map_rows

花费了大约 1.5 小时。

显然是这样的

df = df.with_columns(
        pl.col("pickup_coordinates").list.first().cast(pl.Float32).alias('pickup_longitude'),
        pl.col("pickup_coordinates").list.last().cast(pl.Float32).alias('pickup_latitude'),
        pl.col("dropoff_coordinates").list.first().cast(pl.Float32).alias('dropoff_longitude'),
        pl.col("dropoff_coordinates").list.last().cast(pl.Float32).alias('dropoff_latitude')
    ).with_columns(
        coords = geodesic( (pl.col("pickup_latitude"), pl.col('pickup_longitude')),  (pl.col("dropoff_latitude"), pl.col('dropoff_longitude'))).kilometers
    )

不起作用，因为极坐标尝试对

(pl.col("pickup_latitude"), pl.col('pickup_longitude')

应用逻辑运算

因此，我想了解

map_rows

/

map_elements

是否是我唯一的解决方案，或者是否有不同的解决方法可以加快计算速度。

Answer 1

如https://stackoverflow.com/a/76265233/的答案您可以尝试使用Polars表达式复制

geodesic()

。

另一个潜在的选择可能是 DuckDB，它可以轻松输入/输出 Polars DataFrames。

DuckDB 有一个 SPATIAL 扩展：https://duckdb.org/2023/04/28/spatial.html

duckdb.sql("install spatial") # needed once

如果我增加你的例子进行简单比较：

df_big = pl.concat([df] * 10)

使用您的

map_rows

方法：

(df_big
  .select("pickup_coordinates", "dropoff_coordinates")
  .map_rows(compute_centroid_distance_v2)
  .rename({"map": "centroid_distance"})
)

shape: (98_510, 1)
┌───────────────────┐
│ centroid_distance │
│ ---               │
│ f64               │
╞═══════════════════╡
│ 1.50107           │
│ 0.0               │
│ 0.0               │
│ 3.18019           │
│ 3.652772          │
│ …                 │
│ 2.376629          │
│ 1.440797          │
│ 4.583181          │
│ 0.53954           │
│ 2.589928          │
└───────────────────┘

经过时间：4.52725秒

使用

duckdb

：

duckdb.sql("load spatial")

duckdb.sql("""
from df_big
select
   st_distance_spheroid(
      st_point(
         pickup_coordinates[2]::float, -- NOTE: array indexing is 1-based
         pickup_coordinates[1]::float
      ),
      st_point(
         dropoff_coordinates[2]::float,
         dropoff_coordinates[1]::float
      )
   ) as geodesic
""") .pl()

shape: (98_510, 1)
┌─────────────┐
│ geodesic    │
│ ---         │
│ f64         │
╞═════════════╡
│ 1501.364    │
│ 0.0         │
│ 0.0         │
│ 3180.189287 │
│ 3652.673199 │
│ …           │
│ 2376.786018 │
│ 1440.740571 │
│ 4583.039701 │
│ 539.144276  │
│ 2590.085087 │
└─────────────┘

经过时间：0.10821秒

我对空间数据了解不多，所以我不完全确定为什么输出存在细微差异。

您似乎还可以使用

st_read()

将数据直接加载到 duckdb 中，而不必先使用 Polars 手动将其切碎。

Answer 2

我已经根据我的数据点和提供的解决方案计算了半正弦距离这里

def compute_haversine_disntance(df:pl.DataFrame, R:np.float64, coordinates:dict) -> pl.DataFrame:
    pl.Config.set_fmt_float("full")
    multiplier:float = np.pi/180
    rad_lat1:pl.Expr = (pl.col(coordinates["pickup_points"]).list.last().cast(pl.Float64) * (multiplier))
    rad_lat2:pl.Expr = (pl.col(coordinates["dropoff_points"]).list.last().cast(pl.Float64) * (multiplier))
    rad_lng1:pl.Expr = (pl.col(coordinates["pickup_points"]).list.first().cast(pl.Float64) * (multiplier))
    rad_lng2:pl.Expr = (pl.col(coordinates["dropoff_points"]).list.first().cast(pl.Float64) * (multiplier))
    haversin:pl.Expr = (
        (rad_lat2 - rad_lat1).truediv(2).sin().pow(2) +
        ((rad_lng1.cos() * rad_lng2.cos()) * (rad_lng2 - rad_lng1).truediv(2).sin().pow(2))
    ).cast(pl.Float64)
    df = df.with_columns(
        (
            2 * R * (haversin.sqrt().arcsin())
        ).cast(pl.Float64).alias("haversine_centroid_distance")
    )
    return df

但是，与此计算器这里相比，我的最终结果有一些差异。尽管我的公式与计算器中使用的公式相同，但结果略有不同。例如第一对点：

• lat1：0.7117528862292272（弧度）
• lat2：0.7115465664360616
• lng1：-1.2907933590722993
• lng2：-1.2909438559692195

根据我的计算，

= 的距离为 1.34，而计算器计算出的距离为 1.501（更接近

geodesic

）