具有 GPS 坐标的光线投射算法答案

【问题标题】：Raycasting algorithm with GPS coordinates具有 GPS 坐标的光线投射算法
【发布时间】：2013-01-03 23:43:48
【问题描述】：

我正在使用 Google 地图制作一个小型应用程序，用于确定输入的地址是否属于预定义的服务区域。

用户输入一个地址，PHP 脚本从Geocoding API 获取纬度/经度，并使用构成该区域顶点的一组坐标应用光线投射（取自生成的 KML 文件地图）。

问题是这样的：它在大部分时间都有效，但是服务区域外部的一些地址错误地报告为合格，而该区域内部的其他一些地址则不合格。起初我认为这是 Google 地图的精度问题，但从地理编码服务中的地址生成的坐标是准确的。大概和公式有关吧。

在这里（它基于我在其他地方找到的代码）：

// $points is an array full of Point objects (the vertices), which contain lat/long variables
// $ctr is simply a counter that we will test to see if it's even/odd
for ($i = 0, $j = sizeof($points) - 1; $i < sizeof($points); $j = $i++) {
    $p1 = $points[$i];
    $p2 = $points[$j];
    // $test_point is the lat/long pair of the user's address
    if ($p1->lat < $test_point->lat && $p2->lat >= $test_point->lat ||  $p2->lat < $test_point->lat && $p1->lat >= $test_point->lat)  {
        if ($p1->long + ($test_point->lat - $p1->lat)/($p2->lat - $p1->lat)*($p2->long - $p1->long) < $test_point->long)  
            $ctr++;
    }
}

这里有什么我遗漏的吗？我尝试自己推导出一个公式，并且我在一定程度上理解了这背后的数学原理，但是可以使用来自谷歌地图的 GPS 坐标吗？

对于错误报告的内容似乎没有真正的模式：我测试了靠近边界的地址或服务区角落的地址，但没有运气。另外值得注意的是，这个服务区只是一个城市中相对较小的区域，与州或全国范围内的区域完全不同。

【问题讨论】：

我发现这个算法在这里非常有用：bit.ly/1eJSugg

标签： php gps maps raycasting

【解决方案1】：

嗯....您的第二个 if() 不能补偿任何减法都可能导致负数的事实；仅当坐标严格排序时才有效。

更新：在http://rosettacode.org/wiki/Ray-casting_algorithmn 上有一大堆各种语言的算法详细描述了这个过程（不幸的是，没有 PHP 版本）。您的解决方案中似乎缺少的是选择一个保证在多边形之外的点；因为您正在处理应该很容易的经度/纬度。其次，确保您的多边形是闭合的（即，如果 Google 地图尚未这样做，则从最后一点返回到第一个点）

【讨论】：

纬度/经度可以有负值
这不是重点；该公式在不知道哪个是正确顺序的情况下减去各种变量。在处理顶点和向量时，您必须正确排序，否则您将翻转“内部”或“外部”
好的...但我不一定能取每个减法结果的绝对值并期望得到正确的答案。我怎么知道什么时候是“正确”的顺序？

【解决方案2】：

假设$points 数组包含按顺时针（或逆时针顺序）描述覆盖区域的多边形的角，您的代码对我来说看起来是正确的。基本上，它是在计算与从给定点到第 180 条子午线正东绘制的线相交的多边形边数。

为了清楚起见，我可能会这样重写它：

$p0 = end($points);
foreach ( $points as $p1 ) {
    // ignore edges of constant latitude (yes, this is correct!)
    if ( $p0->lat != $p1->lat ) {
        // scale latitude of $test_point so that $p0 maps to 0 and $p1 to 1:
        $interp = ($test_point->lat - $p0->lat) / ($p1->lat - $p0->lat);
        // does the edge intersect the latitude of $test_point?
        // (note: use >= and < to avoid double-counting exact endpoint hits)
        if ( $interp >= 0 && $interp < 1 ) {
            // longitude of the edge at the latitude of the test point:
            // (could use fancy spherical interpolation here, but for small
            // regions linear interpolation should be fine)
            $long = $interp * $p1->long + (1 - $interp) * $p0->long;
            // is the intersection east of the test point?
            if ( $long < $test_point->long ) {
                // if so, count it:
                $ctr++;
            }
        }
    }
    $p0 = $p1;
}

请注意，如果区域边界越过 180 度子午线，此代码将以各种有趣的方式中断，因此如果您在太平洋中部有任何服务区域，请不要使用它。

如果您仍然遇到问题，请尝试在地图上绘制由$points 数组描述的多边形；您可能会发现它看起来不像您想象的那样，例如如果某些点以错误的顺序列出。

【讨论】：

【解决方案3】：

当光线与形状相切时，此算法存在错误。当它可能发生时，只需在测试点的纬度上添加一个 epsilon（Ilmari 代码的第 3 行）：

if ($test_point->lat == $p0->lat)
    $test_point->lat += 0.0000000001;

另见http://rosettacode.org/wiki/Ray-casting_algorithm（更正网址）。

谢谢。

【讨论】：