如何提高计算hasrsine公式的查询性能？

Question

给出一个有纬度和经度的位置表,哪些位置最接近给定位置？

当然,在地球表面找到距离意味着使用大圆距离,用Haversine公式计算,也称为球面余弦定律公式.

我有以下代码:

 SELECT zip, latitude, longitude, distance
 FROM (
  SELECT z.zip,
         z.latitude, z.longitude,
         p.radius,
         p.distance_unit
             * DEGREES(ACOS(COS(RADIANS(p.latpoint))
             * COS(RADIANS(z.latitude))
             * COS(RADIANS(p.longpoint - z.longitude))
             + SIN(RADIANS(p.latpoint))
             * SIN(RADIANS(z.latitude)))) AS distance
  FROM zip AS z
   JOIN (   /* these are the query parameters */
    SELECT  42.81 AS latpoint, -70.81 AS longpoint,
            50.0  AS radius,  111.045 AS distance_unit
        ) AS p ON 1=1
  WHERE z.latitude
    BETWEEN p.latpoint  - (p.radius / p.distance_unit)
        AND p.latpoint  + (p.radius / p.distance_unit)
   AND z.longitude
    BETWEEN p.longpoint - (p.radius / (p.distance_unit * COS(RADIANS(p.latpoint))))
        AND p.longpoint + (p.radius / (p.distance_unit * COS(RADIANS(p.latpoint))))
) AS d
WHERE distance <= radius

有没有办法提高此查询的性能？

是否有必要使用PostGIS来改进它,或者它只是我的半正式公式的包装？

Answer 1

这个查询永远不会特别快。但是，有一些方法可以对其进行改进。

\n\n

第一：这里不需要半正矢公式。仅当地球曲率是一个重要因素或非常接近两极时，才需要进行校正。但这里的情况并非如此——需要精确计算的最大距离是 12 英里，几乎没有超出地平线。在这个比例上，地球实际上是平坦的，因此毕达哥拉斯定理足以计算距离。

\n\n

1 度纬度约为 69 英里，在 52\xc2\xb0N（荷兰附近）处，1 度经度为 cos(52\xc2\xb0) x 69 = 42.5 英里，因此公式变为：

\n\n

\n

sqrt(pow(69*(lat - $纬度), 2) + pow(42.5*(lng - $经度), 2))

\n

\n\n

第二：我们可以对纬度和经度使用“剪刀测试”。如果某个点在任何基本方向上距目标点超过 12 英里，则它肯定不能位于该点的 12 英里圆圈内。我们可以利用这一事实对纬度和经度进行快速比较，完全跳过距离计算。使用我们上面导出的一度纬度/经度的数字，我们有：

\n\n

\n

WHERE (lat BETWEEN ($latitude - 12/69.0) AND ($latitude + 12/69.0)) \n AND (lng BETWEEN ($longitude - 12/42.5) AND ($longitude + 12/42.5))

\n

\n\n

请注意，这并不能取代全程检查！这只是快速抛出不可能位于正确半径内的点的第一步。通过在纬度或经度上建立索引，数据库服务器可以避免检查数据库中的许多行。

\n

Answer 2

我想规划者会自己重写这个查询，但值得一试。至少它更整洁。

select zip, latitude, longitude, distance
from (
    select z.zip,
           z.latitude, z.longitude,
           p.radius,
           p.distance_unit
               * p.degrees_acos_cos_radians_latpoint
               * cos(radians(z.latitude))
               * cos(radians(p.longpoint - z.longitude))
               + p.sin_radians_latpoint
               * sin(radians(z.latitude)))) as distance
    from
        zip z
        cross join (
            select
                latpoint, longpoint, radius, distance_unit,
                latpoint - radius / distance_unit as lat0,
                latpoint + radius / distance_unit as lat1,
                longpoint - radius / distance_unit * cos(radians(latpoint)) as long0,
                longpoint + radius / distance_unit * cos(radians(latpoint)) as long1,
                sin(radians(latpoint)) as sin_radians_latpoint,
                degrees(acos(cos(radians(latpoint)) as degrees_acos_cos_radians_latpoint
            from (
                values (42.81, -70.81, 50.0, 111.045)
            ) v (latpoint, longpoint, radius, distance_unit)
        ) p
    where
        z.latitude between lat0 and lat1
        and
        z.longitude between long0 and long1
) d
where distance <= radius