使用GPS坐标的射线投射算法

Question

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

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

问题在于：它大多数时候都有效，但有些地址在服务区域之外错误地报告为合格，而其他一些地址则在区域内却不合格。起初我以为这是Google地图的精度问题，但从Geocoding服务中生成的坐标非常准确。这可能与公式有关。

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

// $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坐标是否可行呢？

似乎并没有真正准确的错误报告模式：我测试了靠近边界或服务区角落的地址等情况，但都没有发现规律。值得注意的是，此服务区仅是城市中相对较小的区域，与州或全国范围的区域完全不同。

- John V

---

我发现这里的算法非常有用：http://bit.ly/1eJSugg - user591272

Answer 1

好的...你的第二个if()语句没有考虑到任何减法可能导致负数的情况；它只适用于坐标严格排序的情况。

更新：在http://rosettacode.org/wiki/Ray-casting_algorithmn上，有许多不同语言的算法详细描述了该过程（不幸的是，PHP版本缺失）。您解决方案似乎缺少选择一个保证在多边形外部的点；由于您正在处理经度/纬度，因此应该很容易。其次，请确保您的多边形是封闭的（即从最后一个点返回到第一个点，如果Google Maps尚未执行此操作）。

- JvO

---

经度/纬度可以有负值。 - CSᵠ

这不是重点；该公式在不知道正确顺序的情况下减去各种变量。处理顶点和向量时，必须正确排序，否则会翻转“内部”或“外部”的内容。 - JvO

好的...但我不能简单地对每个减法结果取绝对值并期望得到正确答案。我怎样才能知道哪一个是“正确”的顺序？ - John V

Answer 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数组描述的多边形；您可能会发现它看起来与您想象的不同，例如，如果某些点按错误顺序列出。

- Ilmari Karonen

Answer 3

这个算法存在一个问题，当光线与形状相切时。只需在可能发生的测试点纬度上添加一个 epsilon（Ilmari 代码的第三行）即可解决：

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

也可以查看http://rosettacode.org/wiki/Ray-casting_algorithm（更正的网址）。

谢谢。

- Fil