Haversine 公式单元错误 - PL/SQL

Posted 2023-03-05

【中文标题】Haversine 公式单元错误 - PL/SQL

【英文标题】：Haversine Formula unit errors - PL/SQL

【发布时间】：2017-01-18 17:33:34

【问题描述】：

总体目标：查找指定邮编 30 英里范围内的所有邮编。

数据源：邮政编码列表，经纬度分别为“x”和“y”。

例子：

create or replace view zips as
select 37171 zip, 36.362704 y, -87.30434 x,'Southside' City, 'TN' State from dual 
union
select 37212, 36.133012, -86.802764, 'Nashville', 'TN' from dual 
union
select 37027, 36.00245, -86.791159, 'Brentwood', 'TN' from dual 
union
select 37191, 36.501196, -87.539057, 'Woodlawn', 'TN' from dual 
union
select 37067, 35.928406, -86.805538, 'Franklin', 'TN' from dual ;

我的尝试：我发现了一些看起来像是要去solve all my problems的东西：

CREATE OR REPLACE FUNCTION distance (Lat1 IN NUMBER,
                                     Lon1 IN NUMBER,
                                     Lat2 IN NUMBER,
                                     Lon2 IN NUMBER,
                                     Radius IN NUMBER DEFAULT 3963) RETURN NUMBER IS
 -- Convert degrees to radians
 DegToRad NUMBER := 57.29577951;

BEGIN
  RETURN(NVL(Radius,0) * ACOS((sin(NVL(Lat1,0) / DegToRad) * SIN(NVL(Lat2,0) / DegToRad)) +
        (COS(NVL(Lat1,0) / DegToRad) * COS(NVL(Lat2,0) / DegToRad) *
         COS(NVL(Lon2,0) / DegToRad - NVL(Lon1,0)/ DegToRad))));
END;

但是，它给了我一些奇怪的数字 - 始终太小（我知道从富兰克林到纳什维尔不是 1 英里），但不是我可以确定的恒定因素。 （我用来测试的代码如下 - 我只是将谷歌地图上的距离与这些距离进行了比较）

select b.zip
  ,b.city
  ,b.state
  ,distance(a.x,a.y,b.x,b.y) distance
from zips a, zips b
where a.zip=37067
order by distance;

所以，我认为地理数据组可能有不同的记录纬度和经度的方式，我找到了 Wikipedia's Haversine formula 并将其变成了一个函数

CREATE OR REPLACE FUNCTION distance (Lat1_d IN NUMBER,
                                     Lon1_d IN NUMBER,
                                     Lat2_d IN NUMBER,
                                     Lon2_d IN NUMBER,
                                     Radius IN NUMBER DEFAULT 3959) RETURN NUMBER IS
 -- Convert degrees to radians
 DegToRad NUMBER := .0174532925;
 Lat1 NUMBER := Lat1_d * DegToRad;
 Lon1 NUMBER := Lon1_d * DegToRad;
 Lat2 NUMBER := Lat2_d * DegToRad;
 Lon2 NUMBER := Lon2_d * DegToRad;

BEGIN

  RETURN 2*Radius * 
    asin(
      sqrt(
        power(sin((Lat2-Lat1)/2),2) + cos(Lat1) * cos(Lat2) * power(sin((Lon2-Lon1)/2),2)
      )
    )
        ;
END;

同样的问题 - 结果几乎完全相同，这让我想知道问题出在哪里。我没有做什么转变？我错过了什么？

请注意，我很乐意假设一个平坦的地球并使用毕达哥拉斯距离，因为我只会查看 30-50 英里的每组邮政编码，但我需要弄清楚什么关于我的纬度/经度转换在那里也出错了。

【问题讨论】：

恕我直言，这是一篇很棒的第一篇文章。正确的标签、数据示例、代码和数据格式、参考资料……欢迎来到 SO！

【参考方案1】：

看起来像参数问题。标题是：

CREATE OR REPLACE FUNCTION distance (Lat1 IN NUMBER,
                                     Lon1 IN NUMBER,
                                     Lat2 IN NUMBER,
                                     Lon2 IN NUMBER,

你这样称呼它：

distance(a.x,a.y,b.x,b.y)

但你的观点是：

create or replace view zips as
select 37171 zip, 36.362704 y, -87.30434 x,'Southside' City, 'TN' State from dual 
union
select 37212, 36.133012, -86.802764, 'Nashville', 'TN' from dual 
union
select 37027, 36.00245, -86.791159, 'Brentwood', 'TN' from dual 
union
select 37191, 36.501196, -87.539057, 'Woodlawn', 'TN' from dual 
union
select 37067, 35.928406, -86.805538, 'Franklin', 'TN' from dual ;

根据视图，x 是经度，y 是纬度。所以函数调用应该是：

distance(a.y,a.x,b.y,b.x)

或者，您可以在视图中交换“x”和“y”别名，并保留其他所有内容。

【讨论】：

我听说时不时地用手掌拍打额头很好。动摇事情。感谢您给我机会这样做。

@Maggie，没问题。我也有同样的感觉，并且可能由于多年忘记分号等而在我的头骨上留下了凹痕。

关于体面名字价值的重要一课。如果视图列是 longitude 和 latitude 而不是 x 和 y，您可能会立即发现这一点。