Algorithme de point dans un polygone pour plusieurs polygones

J'ai une carte Google avec un tas de polygones dessus.

Voici un problème qui m'intéresse: Étant donné un point lat, lng, quelle est la meilleure façon de déterminer tous les polygones dans lesquels ce point se trouve?

La manière la plus évidente est d'exécuter un algorithme "point dans le polygone" de manière itérative pour chaque polygone, mais je me demandais s'il existait un algorithme efficace pour répondre à de telles requêtes, surtout si vous avez des milliers de polygones.

Comme pour presque toutes ces questions, l'approche optimale dépend des «cas d'utilisation» et de la façon dont les fonctionnalités sont représentées. Les cas d'utilisation se distinguent généralement par (a) s'il y a plusieurs ou peu d'objets dans chaque couche et (b) si l'une ou l'autre (ou les deux) couches permettent de précalculer certaines structures de données; c'est-à-dire si l'un ou les deux sont suffisamment statiques et immuables pour que l'investissement dans le précalcul en vaille la peine.

Dans le cas présent, cela donne les scénarios suivants. Normalement, les points sont dynamiques: c'est-à-dire qu'ils ne sont pas donnés à l'avance. (S'ils sont disponibles à l'avance ou dans de très grands groupes, certaines optimisations basées sur leur tri seront disponibles.) Soit Q le nombre de points de requête et P le nombre de sommets de polygone .

Données de polygone vectoriel

(1) Peu de points, peu de sommets de polygones au total . Utilisez une procédure de force brute, telle que l' algorithme classique de coup de ligne . Pour toute méthode décente, le coût est O (P * Q), car il faut O (1) de temps pour comparer un point à un bord de polygone et toutes ces comparaisons doivent être faites.

(2) Peut-être plusieurs sommets de polygones, mais ils sont dynamiques: chaque fois qu'un point est utilisé dans la requête, les polygones peuvent tous avoir changé. Utilisez à nouveau un algorithme de force brute. Le coût est toujours O (P * Q), qui sera grand parce que P sera grand, mais cela ne sert à rien. Si les modifications sont petites ou contrôlées ( par exemple , les polygones changent légèrement de forme ou se déplacent simplement lentement), vous pourrez peut-être utiliser une version de la solution suivante et trouver un moyen efficace de mettre à jour les structures de données à mesure que les polygones changent. Ce serait probablement une question de recherche originale.

(3) De nombreux sommets de polygones et polygones statiques (c'est-à-dire que la couche de polygones change rarement). Précalculer une structure de données pour soutenir la recherche (qui pourrait être basé sur un balayage de ligne ou un quadtree algorithme). Le coût du précalcul pour ces algorithmes est O (P * log (P)), mais le coût des requêtes devient O (Q * log (P)), donc le coût total est O ((P + Q) * log ( P)).

Certaines améliorations sont disponibles dans des cas spéciaux , tels que

(a) Tous les polygones sont convexes (le prétraitement des polygones peut être fait plus rapidement ),

(b) Tous les intérieurs de polygone sont disjoints , auquel cas vous pouvez considérer leur union comme étant un seul polygone (ce qui permet des algorithmes efficaces simples, tels que ceux basés sur la triangulation, et

(c) La plupart des polygones ne sont pas très tortueux - c'est-à-dire qu'ils occupent de grandes parties de leurs boîtes englobantes - auquel cas vous pouvez effectuer un test initial basé uniquement sur les boîtes englobantes, puis affiner cette solution. Il s'agit d'une optimisation populaire.

(d) Le nombre de points est important. Les trier pourrait améliorer le timing. Par exemple, lors de la mise en œuvre d'un algorithme de point de polygone de balayage de ligne de gauche à droite, vous devez trier les points sur leur première coordonnée, ce qui vous permet de balayer les points en même temps que vous balayez les bords du polygone. Je ne suis pas au courant qu'une telle optimisation a été publiée. Cependant, celui qui a été publié consiste à effectuer une triangulation contrainte de l'union de tous les points et sommets du polygone: une fois cette triangulation terminée, l'identification des points intérieurs devrait être rapide. Le coût de calcul évoluera comme O (Q * log (Q) + (P + Q) * log (P + Q)).

Données de polygone raster

C'est incroyablement simple: visualisez la couche de polygones comme un raster indicateur binaire (1 = à l'intérieur d'un polygone, 0 = à l'extérieur). (Cela pourrait nécessiter une table de recherche pour convertir les valeurs du raster en indicateurs internes / externes.) Chaque sonde ponctuelle nécessite désormais O (1) pour indexer la cellule raster et lire sa valeur. L'effort total est O (Q).

En général

Une belle solution hybridedans le cas de nombreux polygones vectoriels statiques (cas vectoriel 3 ci-dessus), il s'agit initialement de pixelliser les polygones, peut-être même avec une résolution grossière, en distinguant cette fois toutes les cellules coupant n'importe quelle partie d'une frontière de polygone (donnez-leur une valeur de 2, par exemple) . L'utilisation d'une sonde raster (coût: O (1)) donne généralement une réponse définitive (le point est connu pour être à l'intérieur ou à l'extérieur), mais entraîne parfois une réponse indéfinie (le point tombe dans une cellule à travers laquelle au moins un bord passe), auquel cas la requête vectorielle O (log (P)) la plus coûteuse est effectuée. Cette méthode entraîne des coûts de stockage supplémentaires pour le raster, mais dans de nombreux cas, même un petit raster (un Mo permettra un raster 2000 par 2000 qui stocke les valeurs {0,1,2, null}) peut conférer d'énormes avantages en temps de calcul . Asymptotiquement,

Si vous aviez les boîtes de délimitation de polygone stockées dans quelque chose comme un arbre quadruple, vous pouvez l'utiliser pour déterminer rapidement les polygones à vérifier. Au minimum, vous pouvez simplement voir si le point se trouve à l'intérieur de chaque zone de délimitation de polygone, au lieu de faire un point complet de polygone pour chaque polygone. Personnellement, je configurerais un service Web qui mettrait en cache les polygones en mémoire et utiliserait quelque chose comme JTS ou la suite NetTopology pour faire la requête d'intersection pour moi.

Dans postgis, ST_Intersects utilise des index pour rechercher d'abord si le point se trouve à l'intérieur du cadre de délimitation du polygone, puis vérifie à nouveau s'il se trouve réellement à l'intérieur du polygone. C'est rapide, souvent très rapide.

Si vous avez stocké vos données dans PostGIS, il ne fait aucun doute que la base de données est le bon endroit pour effectuer le calcul. Dans d'autres cas, vous devrez envoyer vos polygones à un programme intermédiaire ou client. Cela, en soi, prendra beaucoup plus de temps que de faire les calculs et d'obtenir simplement les polygones pertinents.

/ Nicklas