Détecter si une coordonnée GPS tombe dans un polygone sur une carte

voici un programme java qui utilise une fonction qui retournera true si une latitude/longitude se trouve à l'intérieur d'un polygone défini par une liste de lat/longs, avec démonstration pour l'état de Floride.

Je ne suis pas sûr qu'il traite du fait que le système lat/long GPS n'est pas un avion de coordonnées x/Y. Pour mes utilisations, j'ai montré que cela fonctionne (je pense que si vous spécifiez suffisamment de points dans la zone de délimitation, il emporte l'effet que la terre est une sphère, et que les lignes droites entre deux points sur la terre ne sont pas une ligne droite de flèche.

tout d'abord spécifier les points qui composent les points d'angle du polygone, il peut avoir des coins concaves et convexes. Les coordonnées que j'utilise ci-dessous retrace le périmètre de l'état de Floride.

méthode coordinate_is_inside_polygon utilise un algorithme que je ne comprends pas très bien. Voici une explication officielle de la source où je l'ai eu:

"... solution transmise par Philippe Reverdy est de calculer la somme des angles fait entre le point d'essai et chaque paire de points constituant le polygone. Si cette somme est 2pi alors le point est un point intérieur, si 0 alors le point est un point extérieur. Cela fonctionne également pour les polygones avec des trous étant donné que le polygone est défini avec un chemin composé d'arêtes coïncidentes entrant et sortant du trou comme c'est la pratique courante dans de nombreux paquets CAO. "

les tests de mon unité montrent qu'il fonctionne de manière fiable, même lorsque la boîte est en forme de " C " ou même la forme d'un Tore. (Mon unité Teste de nombreux points à L'intérieur de la Floride et s'assurer que la fonction retourne vrai. Et je choisis un certain nombre de coordonnées partout ailleurs dans le monde et m'assure qu'il retourne false. Je choisis des endroits dans le monde entier qui pourraient le confondre.

Je ne suis pas sûr que cela fonctionne si la boîte de limite du polygone traverse l'Équateur, Le Méridien Premier, ou n'importe quelle zone où les coordonnées changent de -180 -> 180, -90 -> 90. Ou votre polygone s'enroule autour de la la terre autour de l'axe nord/sud pôles. Pour moi, je n'en ai besoin que pour le périmètre de la Floride. Si vous devez définir un polygone qui enjambe la terre ou qui croise ces lignes, vous pouvez travailler autour de lui en faisant deux polygones, l'un représentant la zone d'un côté du méridien et l'autre représentant la zone de l'autre côté et de tester si votre point est dans l'un de ces points.

voici où j'ai trouvé cet algorithme:Déterminer si un point se trouve à l'intérieur d'une polygon-Solution 2

Exécuter vous-même pour vérifier.

mettez ceci dans un fichier appelé 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's unused, to reject invalid lat/longs.
    if (latitude > -90 && latitude < 90 && 
            longitude > -180 && longitude < 180)
    {
        return true;
    }
    return false;
}
}

la magie démoniaque doit être testée à l'unité. Mettez ça dans un fichier appelé MainTest.java pour vérifier que cela fonctionne pour vous

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.

j'ai utilisé l'IDE eclipse pour faire tourner java en utilisant java 1.6.0. Pour moi tous les tests unitaires passent. Vous devez inclure le fichier JAR junit 4 dans votre classpath ou importer dans Eclipse.

j'ai pensé de la même façon que shab d'abord (sa proposition s'appelle Algorithme De Moulage Des Rayons), mais a des doutes comme Spacedman:

...mais toute la géométrie devra être refait en coordonnées sphériques...

j'ai mis en oeuvre et testé la méthode mathématique correcte pour le faire, E. I. en croisant de grands cercles et en déterminant si l'un des deux points d'intersection se trouve sur les deux arcs. (Note: j'ai suivi les étapes décrites ici, mais j'ai trouvé plusieurs erreurs:sign la fonction est manquante à la fin de l'étape 6 (juste avant arcsin), et le test final est des déchets numériques (comme soustraction est mal conditionnée); utilisez plutôt L_1T >= max(L_1a, L_1b) pour vérifier si S1 est sur le premier arc etc.)

C'est très lent et d'un nombre cauchemar (évalue ~100 fonctions trigonométriques, entre autres); il s'est avéré qu'il n'était pas utilisable dans notre système embarqué système.

There is a trick, though: si la zone considérée est assez petite, il suffit de faire une projection cartographique standard, par exemple projection sphérique de Mercator, de chaque point:

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

Ensuite, vous pouvez appliquer ray-casting, où l'intersection des arcs est vérifiée par cette fonction:

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;
}

si vous avez des coordonnées WGS84 sur la sphère, alors votre polygone divise la sphère en deux zones - comment savoir quelle zone est "à l'intérieur" et laquelle est "à l'extérieur" du polygone? La question est essentiellement dénuée de sens!

par exemple, supposons que le polygone forme la ligne de l'Équateur - l'hémisphère nord est-il " en " ou "en"?

de mémoire, la façon de déterminer si un point se trouve à l'intérieur d'un polygone est d'imaginer tracer une ligne de la position à un point éloigné. Vous comptez ensuite le nombre d'intersections entre la ligne et les segments de ligne du polygone. Si ce nombre est pair, alors il ne se trouve pas dans le polygone. Si elle est fausse, alors il ne situées à l'intérieur du polygone.

en supposant que vous vous occupiez du cas de l'enroulement autour du méridien et de la traversée de l'Équateur (en ajoutant des offsets) - ne pouvez-vous pas traiter cela comme un simple point 2d en polygone ?

Voici l'algorithme écrit en Aller: Il prend les coordonnées de point dans le format [lat, long] et le polygone dans le format [[lat,long], [lat,long]...]. Algorithme va rejoindre le premier et le dernier point du polygone tranche

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
}

Coureur.Code Java dans VB.NET

pour le bénéfice des gens de .NET le même code est mis dans VB.NET. Ont essayé et est assez rapide. Essayé avec 350000 disques, il se termine en quelques minutes. Mais comme l'a dit l'auteur, Je ne suis pas encore à tester des scénarios croisant l'Équateur, multizones, etc.

Utilisation

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

Fonctions

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