Обнаружение того, попадает ли координата GPS в многоугольник на карте

Ответ 1

Вот Java-программа, которая использует функцию, которая будет возвращать true, если широта/долгота найдены внутри многоугольника, определенного списком широты/долготы, с демонстрацией состояния Флориды.

Я не уверен, имеет ли это дело с тем фактом, что система GPS широта/долгота не является координатной плоскостью x/y. Для моего использования я продемонстрировал, что это работает (я думаю, что если вы укажете достаточное количество точек в ограничительной рамке, это размывает эффект того, что земля является сферой, а прямые линии между двумя точками на земле не являются прямой линией стрелки).,

Сначала укажите точки, которые составляют угловые точки многоугольника, он может иметь вогнутые и выпуклые углы. Координаты, которые я использую ниже, отслеживают периметр штата Флорида.

Метод coordinate_is_inside_polygon использует алгоритм, который я не совсем понимаю. Вот официальное объяснение от источника, где я получил это:

"... решение, предложенное Филиппом Реверди, состоит в том, чтобы вычислить сумму углов, образованных между контрольной точкой и каждой парой точек, составляющих многоугольник. Если эта сумма равна 2pi, то точка является внутренней точкой, если 0, то точка является внешней точкой. Это также работает для многоугольников с отверстиями, если полигон определен путем, состоящим из совпадающих ребер в и из отверстия, как это принято во многих пакетах САПР ".

Мои модульные тесты показывают, что он работает надежно, даже когда ограничивающая рамка имеет форму буквы "С" или даже имеет форму тора. (Мои модульные тесты проверяют множество точек во Флориде и проверяют, возвращает ли функция значение true. И я выбираю несколько координат в любой точке мира и проверяю, что она возвращает false. Я выбираю места по всему миру, которые могут его запутать.

Я не уверен, что это сработает, если ограничивающий прямоугольник пересекает экватор, меридиан или любую область, где координаты меняются от -180 → 180, -90 → 90. Или ваш многоугольник окружает землю вокруг Северный/Южный полюса. Для меня это нужно только для работы по периметру Флориды. Если вам нужно определить многоугольник, который охватывает землю или пересекает эти линии, вы можете обойти его, создав два многоугольника, один из которых представляет область на одной стороне меридиана, а другой - область на другой стороне и проверит, является ли ваша точка находится в любой из этих точек.

Вот где я нашел этот алгоритм: Определение, лежит ли точка внутри многоугольника - Решение 2

Запустите его для себя, чтобы дважды проверить это.

Поместите это в файл с именем Runner.java

import java.util.ArrayList;
public class Runner
{
    public static double PI = 3.14159265;
    public static double TWOPI = 2*PI;
    public static void main(String[] args) {
    ArrayList<Double> lat_array = new ArrayList<Double>();
    ArrayList<Double> long_array = new ArrayList<Double>();

    //This is the polygon bounding box, if you plot it, 
    //you'll notice it is a rough tracing of the parameter of 
    //the state of Florida starting at the upper left, moving 
    //clockwise, and finishing at the upper left corner of florida.

    ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();
    polygon_lat_long_pairs.add("31.000213,-87.584839");  
    //lat/long of upper left tip of florida.
    polygon_lat_long_pairs.add("31.009629,-85.003052");
    polygon_lat_long_pairs.add("30.726726,-84.838257");
    polygon_lat_long_pairs.add("30.584962,-82.168579");
    polygon_lat_long_pairs.add("30.73617,-81.476441");  
    //lat/long of upper right tip of florida.
    polygon_lat_long_pairs.add("29.002375,-80.795288");
    polygon_lat_long_pairs.add("26.896598,-79.938355");
    polygon_lat_long_pairs.add("25.813738,-80.059204");
    polygon_lat_long_pairs.add("24.93028,-80.454712");
    polygon_lat_long_pairs.add("24.401135,-81.817017");
    polygon_lat_long_pairs.add("24.700927,-81.959839");
    polygon_lat_long_pairs.add("24.950203,-81.124878");
    polygon_lat_long_pairs.add("26.0015,-82.014771");
    polygon_lat_long_pairs.add("27.833247,-83.014527");
    polygon_lat_long_pairs.add("28.8389,-82.871704");
    polygon_lat_long_pairs.add("29.987293,-84.091187");
    polygon_lat_long_pairs.add("29.539053,-85.134888");
    polygon_lat_long_pairs.add("30.272352,-86.47522");
    polygon_lat_long_pairs.add("30.281839,-87.628784");

    //Convert the strings to doubles.       
    for(String s : polygon_lat_long_pairs){
        lat_array.add(Double.parseDouble(s.split(",")[0]));
        long_array.add(Double.parseDouble(s.split(",")[1]));
    }

   //prints TRUE true because the lat/long passed in is
    //inside the bounding box.
    System.out.println(coordinate_is_inside_polygon(
            25.7814014D,-80.186969D,
            lat_array, long_array));

    //prints FALSE because the lat/long passed in 
    //is Not inside the bounding box.
    System.out.println(coordinate_is_inside_polygon(
            25.831538D,-1.069338D, 
            lat_array, long_array));

}
public static boolean coordinate_is_inside_polygon(
    double latitude, double longitude, 
    ArrayList<Double> lat_array, ArrayList<Double> long_array)
{       
       int i;
       double angle=0;
       double point1_lat;
       double point1_long;
       double point2_lat;
       double point2_long;
       int n = lat_array.size();

       for (i=0;i<n;i++) {
          point1_lat = lat_array.get(i) - latitude;
          point1_long = long_array.get(i) - longitude;
          point2_lat = lat_array.get((i+1)%n) - latitude; 
          //you should have paid more attention in high school geometry.
          point2_long = long_array.get((i+1)%n) - longitude;
          angle += Angle2D(point1_lat,point1_long,point2_lat,point2_long);
       }

       if (Math.abs(angle) < PI)
          return false;
       else
          return true;
}

public static double Angle2D(double y1, double x1, double y2, double x2)
{
   double dtheta,theta1,theta2;

   theta1 = Math.atan2(y1,x1);
   theta2 = Math.atan2(y2,x2);
   dtheta = theta2 - theta1;
   while (dtheta > PI)
      dtheta -= TWOPI;
   while (dtheta < -PI)
      dtheta += TWOPI;

   return(dtheta);
}

public static boolean is_valid_gps_coordinate(double latitude, 
    double longitude)
{
    //This is a bonus function, it unused, to reject invalid lat/longs.
    if (latitude > -90 && latitude < 90 && 
            longitude > -180 && longitude < 180)
    {
        return true;
    }
    return false;
}
}

Магия демона должна быть проверена юнитом. Поместите это в файл с именем MainTest.java, чтобы убедиться, что он работает для вас.

import java.util.ArrayList;
import org.junit.Test;
import static org.junit.Assert.*;

public class MainTest {
@Test
public void test_lat_long_in_bounds(){
    Runner r = new Runner();
    //These make sure the lat/long passed in is a valid gps 
    //lat/long coordinate.  These should be valid. 
    assertTrue(r.is_valid_gps_coordinate(25, -82));
    assertTrue(r.is_valid_gps_coordinate(-25, -82));
    assertTrue(r.is_valid_gps_coordinate(25, 82));
    assertTrue(r.is_valid_gps_coordinate(-25, 82));
    assertTrue(r.is_valid_gps_coordinate(0, 0));
    assertTrue(r.is_valid_gps_coordinate(89, 179));
    assertTrue(r.is_valid_gps_coordinate(-89, -179));
    assertTrue(r.is_valid_gps_coordinate(89.999, 179));
    //If your bounding box crosses the equator or prime meridian, 
    then you have to test for those situations still work.
}
@Test
public void realTest_for_points_inside()
{
    ArrayList<Double> lat_array = new ArrayList<Double>();
    ArrayList<Double> long_array = new ArrayList<Double>();
    ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();
    //upper left tip of florida.
    polygon_lat_long_pairs.add("31.000213,-87.584839");
    polygon_lat_long_pairs.add("31.009629,-85.003052");
    polygon_lat_long_pairs.add("30.726726,-84.838257");
    polygon_lat_long_pairs.add("30.584962,-82.168579");
    polygon_lat_long_pairs.add("30.73617,-81.476441");  
    //upper right tip of florida.
    polygon_lat_long_pairs.add("29.002375,-80.795288");
    polygon_lat_long_pairs.add("26.896598,-79.938355");
    polygon_lat_long_pairs.add("25.813738,-80.059204");
    polygon_lat_long_pairs.add("24.93028,-80.454712");
    polygon_lat_long_pairs.add("24.401135,-81.817017");
    polygon_lat_long_pairs.add("24.700927,-81.959839");
    polygon_lat_long_pairs.add("24.950203,-81.124878");
    polygon_lat_long_pairs.add("26.0015,-82.014771");
    polygon_lat_long_pairs.add("27.833247,-83.014527");
    polygon_lat_long_pairs.add("28.8389,-82.871704");
    polygon_lat_long_pairs.add("29.987293,-84.091187");
    polygon_lat_long_pairs.add("29.539053,-85.134888");
    polygon_lat_long_pairs.add("30.272352,-86.47522");
    polygon_lat_long_pairs.add("30.281839,-87.628784");

    for(String s : polygon_lat_long_pairs){
        lat_array.add(Double.parseDouble(s.split(",")[0]));
        long_array.add(Double.parseDouble(s.split(",")[1]));
    }

    Runner r = new Runner();
    ArrayList<String> pointsInside = new ArrayList<String>();
    pointsInside.add("30.82112,-87.255249");
    pointsInside.add("30.499804,-86.8927");
    pointsInside.add("29.96826,-85.036011");
    pointsInside.add("30.490338,-83.981323");
    pointsInside.add("29.825395,-83.344116");
    pointsInside.add("30.215406,-81.828003");
    pointsInside.add("29.299813,-82.728882");
    pointsInside.add("28.540135,-81.212769");
    pointsInside.add("27.92065,-82.619019");
    pointsInside.add("28.143691,-81.740113");
    pointsInside.add("27.473186,-80.718384");
    pointsInside.add("26.769154,-81.729126");
    pointsInside.add("25.853292,-80.223999");
    pointsInside.add("25.278477,-80.707398");
    pointsInside.add("24.571105,-81.762085");   //bottom tip of keywest
    pointsInside.add("24.900388,-80.663452");
    pointsInside.add("24.680963,-81.366577");

    for(String s : pointsInside)
    {
        assertTrue(r.coordinate_is_inside_polygon(
            Double.parseDouble(s.split(",")[0]), 
            Double.parseDouble(s.split(",")[1]), 
            lat_array, long_array));
    }
}

@Test
public void realTest_for_points_outside()
{
    ArrayList<Double> lat_array = new ArrayList<Double>();
    ArrayList<Double> long_array = new ArrayList<Double>();

    ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();
    //upper left tip, florida.
    polygon_lat_long_pairs.add("31.000213,-87.584839");
    polygon_lat_long_pairs.add("31.009629,-85.003052");
    polygon_lat_long_pairs.add("30.726726,-84.838257");
    polygon_lat_long_pairs.add("30.584962,-82.168579");
    polygon_lat_long_pairs.add("30.73617,-81.476441");
    //upper right tip, florida.
    polygon_lat_long_pairs.add("29.002375,-80.795288");
    polygon_lat_long_pairs.add("26.896598,-79.938355");
    polygon_lat_long_pairs.add("25.813738,-80.059204");
    polygon_lat_long_pairs.add("24.93028,-80.454712");
    polygon_lat_long_pairs.add("24.401135,-81.817017");
    polygon_lat_long_pairs.add("24.700927,-81.959839");
    polygon_lat_long_pairs.add("24.950203,-81.124878");
    polygon_lat_long_pairs.add("26.0015,-82.014771");
    polygon_lat_long_pairs.add("27.833247,-83.014527");
    polygon_lat_long_pairs.add("28.8389,-82.871704");
    polygon_lat_long_pairs.add("29.987293,-84.091187");
    polygon_lat_long_pairs.add("29.539053,-85.134888");
    polygon_lat_long_pairs.add("30.272352,-86.47522");
    polygon_lat_long_pairs.add("30.281839,-87.628784");

    for(String s : polygon_lat_long_pairs)
    {
        lat_array.add(Double.parseDouble(s.split(",")[0]));
        long_array.add(Double.parseDouble(s.split(",")[1]));
    }

    Runner r = new Runner();

    ArrayList<String> pointsOutside = new ArrayList<String>();
    pointsOutside.add("31.451159,-87.958374");
    pointsOutside.add("31.319856,-84.607544");
    pointsOutside.add("30.868282,-84.717407");
    pointsOutside.add("31.338624,-81.685181");
    pointsOutside.add("29.452991,-80.498657");
    pointsOutside.add("26.935783,-79.487915");
    pointsOutside.add("25.159207,-79.916382");
    pointsOutside.add("24.311058,-81.17981");
    pointsOutside.add("25.149263,-81.838989");
    pointsOutside.add("27.726326,-83.695679");
    pointsOutside.add("29.787263,-87.024536");
    pointsOutside.add("29.205877,-62.102052");
    pointsOutside.add("14.025751,-80.690919");
    pointsOutside.add("29.029276,-90.805666");
    pointsOutside.add("-12.606032,-70.151369");
    pointsOutside.add("-56.520716,-172.822269");
    pointsOutside.add("-75.89666,9.082024");
    pointsOutside.add("-24.078567,142.675774");
    pointsOutside.add("84.940737,177.480462");
    pointsOutside.add("47.374545,9.082024");
    pointsOutside.add("25.831538,-1.069338");
    pointsOutside.add("0,0");

    for(String s : pointsOutside){
        assertFalse(r.coordinate_is_inside_polygon(
            Double.parseDouble(s.split(",")[0]),
            Double.parseDouble(s.split(",")[1]), lat_array, long_array));
    }
}
}
//The list of lat/long inside florida bounding box all return true.
//The list of lat/long outside florida bounding box all return false.

Я использовал eclipse IDE, чтобы запустить Java с использованием Java 1.6.0. Для меня все юнит-тесты пройдены. Вам нужно включить jar файл junit 4 в ваш путь к классам или импортировать его в Eclipse.

Ответ 2

Сначала я думал так же, как и Шэб (его предложение называется " Алгоритм приведения к лучу"), но у меня были такие мысли, как у Spacedman:

... но всю геометрию придется переделывать в сферических координатах...

Я реализовал и протестировал математически правильный способ сделать это, пересекая большие круги и определяя, находится ли одна из двух пересекающихся точек на обеих дугах. (Примечание: я выполнил шаги, описанные здесь, но обнаружил несколько ошибок: функция sign отсутствует в конце шага 6 (непосредственно перед arcsin), а последний тест - числовой мусор (так как вычитание плохо обусловлено); используйте скорее L_1T >= max(L_1a, L_1b) чтобы проверить, находится ли S1 на первой дуге и т.д.)

Это также очень медленный и численный кошмар (среди прочего, оценивается ~ 100 тригонометрических функций); он оказался непригодным для использования в наших встроенных системах.

Однако есть хитрость: если рассматриваемая область достаточно мала, просто сделайте стандартную картографическую проекцию, например сферическую проекцию Меркатора, для каждой точки:

// latitude, longitude in radians
x = longitude;
y = log(tan(pi/4 + latitude/2));

Затем вы можете применить приведение лучей, где пересечение дуг проверяется этой функцией:

public bool ArcsIntersecting(double x1, double y1, double x2, double y2, 
  double x3, double y3, double x4, double y4)
    {

    double vx1 = x2 - x1;
    double vy1 = y2 - y1;

    double vx2 = x4 - x3;
    double vy2 = y4 - y3;

    double denom = vx1 * vy2 - vx2 * vy1;

    if (denom == 0) { return false; } // edges are parallel

    double t1 = (vx2 * (y1 - y3) - vy2 * (x1 - x3)) / denom;

    double t2;

    if (vx2 != 0) { t2 = (x1 - x3 + t1 * vx1) / vx2; }
    else if (vy2 != 0) { t2 = (y1 - y3 + t1 * vy1) / vy2; }
    else { return false; } // edges are matching

    return min(t1, t2) >= 0 && max(t1, t2) <= 1;
}

Ответ 3

Если у вас есть координаты WGS84 в сфере, ваш многоугольник делит сферу на две области - как мы узнаем, какая область "внутри" и которая "снаружи" многоугольника? Вопрос практически бессмысленен!

Например, предположим, что многоугольник образовал линию экватора - северное полушарие "в" или "вне"?

Ответ 4

Из памяти способ определить, находится ли точка в многоугольнике, - это представить, как рисовать линию от позиции до некоторой отдаленной точки. Затем вы подсчитываете количество пересечений между линией и отрезками линии многоугольника. Если он считается четным, то он не лежит внутри полигона. Если оно ложно, то оно лежит внутри многоугольника.

Ответ 5

Предполагая, что вы обрабатываете случай обертывания вокруг меридиана и пересечения экватора (путем добавления смещений) - не можете ли вы просто рассматривать это как простую 2d-точку в многоугольнике?

Ответ 6

Вот алгоритм, написанный в Go: Он принимает координаты точки в формате [lat, long] и polygon в формате [[lat, long], [lat, long]...]. Алгоритм присоединяется к первой и последней точке в многоугольном срезе

import "math"

// ContainsLocation determines whether the point is inside the polygon
func ContainsLocation(point []float64, polygon [][]float64, geodesic 
   bool) bool {
    size := len(polygon)
    if size == 0 {
        return false
    }

    var (
        lat2, lng2, dLng3 float64
    )

    lat3 := toRadians(point[0])
    lng3 := toRadians(point[1])
    prev := polygon[size-1]
    lat1 := toRadians(prev[0])
    lng1 := toRadians(prev[1])
    nIntersect := 0

    for _, v := range polygon {
        dLng3 = wrap(lng3-lng1, -math.Pi, math.Pi)
        // Special case: point equal to vertex is inside.
        if lat3 == lat1 && dLng3 == 0 {
            return true
        }
        lat2 = toRadians(v[0])
        lng2 = toRadians(v[1])
        // Offset longitudes by -lng1.
        if intersects(lat1, lat2, wrap(lng2-lng1, -math.Pi, math.Pi), lat3, dLng3, geodesic) {
            nIntersect++
        }
        lat1 = lat2
        lng1 = lng2
    }

    return (nIntersect & 1) != 0
}

func toRadians(p float64) float64 {
    return p * (math.Pi / 180.0)
}

func wrap(n, min, max float64) float64 {
    if n >= min && n < max {
        return n
    }
    return mod(n-min, max-min) + min
}

func mod(x, m float64) float64 {
    return math.Remainder(math.Remainder(x, m)+m, m)
}

func intersects(lat1, lat2, lng2, lat3, lng3 float64, geodesic bool) bool {
    // Both ends on the same side of lng3.
    if (lng3 >= 0 && lng3 >= lng2) || (lng3 < 0 && lng3 < lng2) {
        return false
    }
    // Point is South Pole.
    if lat3 <= -math.Pi/2 {
        return false
    }
    // Any segment end is a pole.
    if lat1 <= -math.Pi/2 || lat2 <= -math.Pi/2 || lat1 >= math.Pi/2 || lat2 >= math.Pi/2 {
        return false
    }
    if lng2 <= -math.Pi {
        return false
    }
    linearLat := (lat1*(lng2-lng3) + lat2*lng3) / lng2
    // Northern hemisphere and point under lat-lng line.
    if lat1 >= 0 && lat2 >= 0 && lat3 < linearLat {
        return false
    }
    // Southern hemisphere and point above lat-lng line.
    if lat1 <= 0 && lat2 <= 0 && lat3 >= linearLat {
        return true
    }
    // North Pole.
    if lat3 >= math.Pi/2 {
        return true
    }

    // Compare lat3 with latitude on the GC/Rhumb segment corresponding to lng3.
    // Compare through a strictly-increasing function (tan() or mercator()) as convenient.
    if geodesic {
        return math.Tan(lat3) >= tanLatGC(lat1, lat2, lng2, lng3)
    }
    return mercator(lat3) >= mercatorLatRhumb(lat1, lat2, lng2, lng3)
}

func tanLatGC(lat1, lat2, lng2, lng3 float64) float64 {
    return (math.Tan(lat1)*math.Sin(lng2-lng3) + math.Tan(lat2)*math.Sin(lng3)) / math.Sin(lng2)
}

func mercator(lat float64) float64 {
    return math.Log(math.Tan(lat*0.5 + math.Pi/4))
}

func mercatorLatRhumb(lat1, lat2, lng2, lng3 float64) float64 {
    return (mercator(lat1)*(lng2-lng3) + mercator(lat2)*lng3) / lng2
}

Ответ 7

Код Runner.Java в VB.NET

В интересах .NET-пользователей тот же код помещается в VB.NET. Попробовали и довольно быстро. Пробовал 350000 записей, он заканчивается всего за несколько минут. Но, как сказал автор, я еще не тестировал сценарии, пересекающие экватор, мультизоны и т.д.

'Использование

If coordinate_is_inside_polygon(CurLat, CurLong, Lat_Array, Long_Array) Then
    MsgBox("Location " & CurLat & "," & CurLong & " is within polygon boundary")
Else
    MsgBox("Location " & CurLat & "," & CurLong & " is NOT within polygon boundary")
End If

'Функции

Public Function coordinate_is_inside_polygon(ByVal latitude As Double, ByVal longitude As Double, ByVal lat_array() As Double, ByVal long_array() As Double) As Boolean
    Dim i As Integer
    Dim angle As Double = 0
    Dim point1_lat As Double
    Dim point1_long As Double
    Dim point2_lat As Double
    Dim point2_long As Double
    Dim n As Integer = lat_array.Length()
    For i = 0 To n - 1
        point1_lat = lat_array(i) - latitude
        point1_long = long_array(i) - longitude
        point2_lat = lat_array((i + 1) Mod n) - latitude
        point2_long = long_array((i + 1) Mod n) - longitude
        angle += Angle2D(point1_lat, point1_long, point2_lat, point2_long)
    Next

    If Math.Abs(angle) < PI Then Return False Else Return True
End Function

Public Function Angle2D(ByVal y1 As Double, ByVal x1 As Double, ByVal y2 As Double, ByVal x2 As Double) As Double
    Dim dtheta, theta1, theta2 As Double
    theta1 = Math.Atan2(y1, x1)
    theta2 = Math.Atan2(y2, x2)
    dtheta = theta2 - theta1
    While dtheta > PI
         dtheta -= TWOPI
    End While

    While dtheta < -PI
        dtheta += TWOPI
    End While
    Return (dtheta)
End Function



 Public Function is_valid_gps_coordinate(ByVal latitude As Double, ByVal longitude As Double) As Boolean
        If latitude > -90 AndAlso latitude < 90 AndAlso longitude > -180 AndAlso longitude < 180 Then
            Return True
        End If

        Return False
End Function

Ответ 8

Версия JavaScript -

{
const PI = 3.14159265;
const TWOPI = 2*PI;
function isCoordinateInsidePitch(latitude, longitude, latArray, longArray)
{       
       let angle=0;
       let p1Lat;
       let p1Long;
       let p2Lat;
       let p2Long;
       let n = latArray.length;

       for (let i = 0; i < n; i++) {
          p1Lat = latArray[i] - latitude;
          p1Long = longArray[i] - longitude;
          p2Lat = latArray[(i+1)%n] - latitude;
          p2Long = longArray[(i+1)%n] - longitude;
          angle += angle2D(p1Lat,p1Long,p2Lat,p2Long);
       }

       return !(Math.abs(angle) < PI);
}

function angle2D(y1, x1, y2, x2)
{
   let dtheta,theta1,theta2;

   theta1 = Math.atan2(y1,x1);
   theta2 = Math.atan2(y2,x2);
   dtheta = theta2 - theta1;
   while (dtheta > PI)
      dtheta -= TWOPI;
   while (dtheta < -PI)
      dtheta += TWOPI;

   return dtheta;
}

function isValidCoordinate(latitude,longitude)
{
    return (
    latitude !== '' && longitude !== '' && !isNaN(latitude) 
     && !isNaN(longitude) && latitude > -90 &&
     latitude < 90 && longitude > -180 && longitude < 180
     )
}
let latArray = [32.10458, 32.10479, 32.1038, 32.10361];
let longArray = [34.86448, 34.86529, 34.86563, 34.86486];
// true
console.log(isCoordinateInsidePitch(32.104447, 34.865108,latArray, longArray));
// false
// isCoordinateInsidePitch(32.104974, 34.864576,latArray, longArray);
// true
// isValidCoordinate(0, 0)
// true
// isValidCoordinate(32.104974, 34.864576)
}