Suite

Algorithme de point dans le 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 se trouve ce point ?

Le moyen évident consiste à exécuter un algorithme "point dans le polygone" de manière itérative pour chaque polygone, mais je me demandais s'il existe 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 manière dont les fonctionnalités sont représentées. Les cas d'utilisation se distinguent généralement par (a) s'il y a beaucoup ou peu d'objets dans chaque couche et (b) si l'une (ou les deux) couches permet le précalcul de 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 en très grands groupes, des optimisations basées sur leur tri seront disponibles.) Soit Q être le nombre de points de requête et P être le nombre de polygone sommets.

Données de polygone vectoriel

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

(2) Peut-être de nombreux 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), ce qui sera important car P sera grand, mais il n'y a pas d'aide à cela. Si les changements sont mineurs ou contrôlés (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. Cela ferait probablement l'objet d'une recherche originale.

(3) De nombreux sommets de polygones et polygones statiques (c'est-à-dire que la couche de polygones changera rarement). Précalculez une structure de données pour prendre en charge la recherche (qui peut être basée sur un balayage de ligne ou un algorithme quadtree). 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 cas spéciaux, tel que

(une) Tous les polygones sont convexes (le prétraitement des polygones peut se faire plus rapidement),

(b) Tous les intérieurs de polygones sont disjoints, auquel cas vous pouvez considérer leur union comme étant un seul polygone (ce qui permet des algorithmes simples et efficaces, 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 une grande partie de leurs cadres englobants -- auquel cas vous pouvez effectuer un test initial basé uniquement sur les cadres englobants, puis affiner cette solution. Il s'agit d'une optimisation populaire.

(ré) Le nombre de points est grand. Les trier pourrait améliorer le timing. Par exemple, lors de la mise en œuvre d'un algorithme de point dans un polygone de balayage de ligne de gauche à droite, vous triez 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 ait été publiée. Une qui a été publiée, cependant, 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 sera égal à O(Q*log(Q) + (P+Q)*log(P+Q)).

Données de polygone raster

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

En général

Une belle solution hybride dans le cas de nombreux polygones vectoriels statiques (cas vectoriel 3 ci-dessus) est d'abord de rastériser 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 limite de polygone (donnez-leur une valeur de 2, disons) . 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 de 2000 par 2000 qui stocke les valeurs {0,1,2,null}) peut conférer d'énormes avantages en temps de calcul . Asymptotiquement, l'effort de calcul est le même que pour une solution vectorielle, mais en pratique il est de O(Q + P*log(P)) et peut-être aussi faible que O(Q+P) (obtenu en utilisant une résolution très fine pour le raster et en utilisant des méthodes de force brute pour les très rares requêtes vectorielles qui doivent être effectuées).


Si vous aviez stocké les cadres de délimitation des polygones dans quelque chose comme un arbre quad, vous pourriez l'utiliser pour déterminer rapidement quels polygones vérifier. Au minimum, vous pouvez simplement voir si le point est à l'intérieur de chaque boîte englobante de polygone au lieu de faire un point complet dans un 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 effectuer la requête d'intersection pour moi.


Dans postgis, ST_Intersects utilise des index pour d'abord trouver si le point est à l'intérieur de la boîte englobante du polygone, puis vérifie à nouveau s'il se trouve vraiment à 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