Comment générer une grille irrégulière contenant n points minimum?

20

Étant donné un large échantillon (~ 1 million) de points inégalement répartis - est-il possible de générer une grille irrégulière (en taille, mais pourrait également être de forme irrégulière si cela est possible?) Qui contiendra la quantité minimale spécifiée de n points?

Cela a moins d'importance pour moi si les «cellules» générées d'une telle grille contiennent exactement n nombre de points ou au moins n points.

Je connais des solutions comme genvecgrid dans ArcGIS ou Create Grid Layer dans QGIS / mmgis, mais elles créeront toutes des grilles régulières qui se traduiront par une sortie avec des cellules vides (problème plus petit - je pourrais simplement les éliminer) ou des cellules avec un nombre de points moins que n (plus gros problème car j'aurais besoin d'une solution pour agréger ces cellules, en utilisant probablement des outils d' ici ?).

Je suis allé sur Google en vain et je suis ouvert aux solutions commerciales (ArcGIS et extensions) ou gratuites (Python, PostGIS, R).

radek
la source
1
Dans quelle mesure le réseau doit-il être "régulier"? Je me demande si vous pouvez faire un regroupement hiérarchique, puis couper le dendrogramme pour répondre à vos besoins (bien que cela s'étende probablement ce qui serait défini comme une configuration spatiale régulière). La documentation de CrimeStat contient de bons exemples de ce type de clustering .
Andy W
5
Pourriez-vous expliquer exactement ce que vous entendez par "grille irrégulière"? Cela ressemble à un oxymore :-). Plus précisément, quel serait le but de cet exercice? Notez également que des critères ou des contraintes supplémentaires sont probablement nécessaires: après tout, si vous dessinez un carré autour des 1 million de points, il pourrait être considéré comme faisant partie d'une grille et il en contiendrait plus de n . Vous ne vous soucieriez probablement pas de cette solution triviale, cependant: mais pourquoi pas, exactement?
whuber
@AndyW Merci. Bonne idée et mérite d'être explorée. Aura un coup d'oeil. Taille et forme de « grille » est d' une importance secondaire pour moi - priorité ( en raison de la confidentialité des données) est de « cacher » n caractéristiques derrière un
radek
@whuber Merci aussi. Je suis d'accord - mais je ne savais pas trop comment nommer un tel partitionnement. Comme mentionné ci-dessus - ma principale motivation est la confidentialité des données. Ayant cinq emplacements de points (que je ne peux pas montrer sur la carte finale), je voudrais les représenter par zone les couvrant; et obtenir moyenne / médiane / etc. valeur pour cela. Je suis d'accord qu'il serait possible de dessiner un rectangle ou une coque convexe les représentant tous - ce serait la protection ultime de la confidentialité des données, je suppose? ;] Cependant - il serait plus utile de le représenter par des formes délimitant, disons 10 entités. Ensuite - je peux toujours conserver le modèle spatial.
radek
1
IMO étant donné votre description, j'utiliserais un certain type d'interpolation et afficherais une carte raster (peut-être qu'une bande passante adaptative de la taille de votre N minimal serait suffisante pour lisser les données). En ce qui concerne CrimeStat, les plus gros fichiers que j'ai utilisés étaient environ 100 000 cas, je crois (et le clustering prendrait certainement du temps). Il est probable que vous puissiez effectuer une pré-généralisation de vos données pour les représenter comme moins de cas et toujours obtenir des résultats souhaitables pour tout ce que vous voulez. C'est un programme vraiment simple, je suggère de prendre quelques minutes pour l'essayer et voir par vous-même.
Andy W

Réponses:

26

Je vois que MerseyViking a recommandé un quadtree . J'allais suggérer la même chose et pour l'expliquer, voici le code et un exemple. Le code est écrit Rmais doit facilement être porté sur, par exemple, Python.

L'idée est remarquablement simple: diviser les points environ en deux dans la direction x, puis diviser récursivement les deux moitiés le long de la direction y, en alternant les directions à chaque niveau, jusqu'à ce que plus aucune division ne soit souhaitée.

Étant donné que l'intention est de masquer les emplacements réels des points, il est utile d'introduire un certain caractère aléatoire dans les divisions . Un moyen simple et rapide de le faire est de diviser à un ensemble de quantiles une petite quantité aléatoire à partir de 50%. De cette façon (a) les valeurs de division sont très peu susceptibles de coïncider avec les coordonnées des données, de sorte que les points tomberont uniquement dans les quadrants créés par le partitionnement, et (b) les coordonnées des points seront impossibles à reconstruire avec précision à partir du quadtree.

Parce que l'intention est de maintenir une quantité minimale kde nœuds dans chaque feuille de quadtree, nous implémentons une forme restreinte de quadtree. Il prendra en charge (1) les points de regroupement en groupes ayant entre ket 2 * k-1 éléments chacun et (2) la cartographie des quadrants.

Ce Rcode crée un arbre de nœuds et de feuilles terminales, les distinguant par classe. L'étiquetage de classe accélère le post-traitement tel que le traçage, illustré ci-dessous. Le code utilise des valeurs numériques pour les identifiants. Cela fonctionne jusqu'à des profondeurs de 52 dans l'arborescence (en utilisant des doubles; si des entiers longs non signés sont utilisés, la profondeur maximale est de 32). Pour les arbres plus profonds (qui sont très peu probables dans n'importe quelle application, car au moins k* 2 ^ 52 points seraient impliqués), les identifiants devraient être des chaînes.

quadtree <- function(xy, k=1) {
  d = dim(xy)[2]
  quad <- function(xy, i, id=1) {
    if (length(xy) < 2*k*d) {
      rv = list(id=id, value=xy)
      class(rv) <- "quadtree.leaf"
    }
    else {
      q0 <- (1 + runif(1,min=-1/2,max=1/2)/dim(xy)[1])/2 # Random quantile near the median
      x0 <- quantile(xy[,i], q0)
      j <- i %% d + 1 # (Works for octrees, too...)
      rv <- list(index=i, threshold=x0, 
                 lower=quad(xy[xy[,i] <= x0, ], j, id*2), 
                 upper=quad(xy[xy[,i] > x0, ], j, id*2+1))
      class(rv) <- "quadtree"
    }
    return(rv)
  }
  quad(xy, 1)
}

Notez que la conception récursive de division et de conquête de cet algorithme (et, par conséquent, de la plupart des algorithmes de post-traitement) signifie que l'exigence de temps est O (m) et que l'utilisation de la RAM est O (n) où mest le nombre de cellules et nest le nombre de points. mest proportionnel à ndivisé par le nombre minimum de points par cellule,k. Ceci est utile pour estimer les temps de calcul. Par exemple, s'il faut 13 secondes pour partitionner n = 10 ^ 6 points en cellules de 50-99 points (k = 50), m = 10 ^ 6/50 = 20000. Si vous voulez plutôt partitionner en 5-9 points par cellule (k = 5), m est 10 fois plus grand, donc le chronométrage monte à environ 130 secondes. (Parce que le processus de division d'un ensemble de coordonnées autour de leur milieu devient plus rapide à mesure que les cellules deviennent plus petites, le temps réel n'était que de 90 secondes.) Pour aller jusqu'à k = 1 point par cellule, cela prendra environ six fois plus longtemps encore, ou neuf minutes, et nous pouvons nous attendre à ce que le code soit un peu plus rapide que cela.

Avant d'aller plus loin, générons des données intéressantes espacées de manière irrégulière et créons leur quadtree restreint (temps écoulé de 0,29 seconde):

Quadtree

Voici le code pour produire ces tracés. Il exploite Rle polymorphisme de: points.quadtreesera appelé chaque fois que la pointsfonction est appliquée à un quadtreeobjet, par exemple. La puissance de ceci est évidente dans l'extrême simplicité de la fonction pour colorer les points selon leur identifiant de cluster:

points.quadtree <- function(q, ...) {
  points(q$lower, ...); points(q$upper, ...)
}
points.quadtree.leaf <- function(q, ...) {
  points(q$value, col=hsv(q$id), ...)
}

Le tracé de la grille elle-même est un peu plus délicat car il nécessite un écrêtage répété des seuils utilisés pour le partitionnement en quadtree, mais la même approche récursive est simple et élégante. Utilisez une variante pour construire des représentations polygonales des quadrants si vous le souhaitez.

lines.quadtree <- function(q, xylim, ...) {
  i <- q$index
  j <- 3 - q$index
  clip <- function(xylim.clip, i, upper) {
    if (upper) xylim.clip[1, i] <- max(q$threshold, xylim.clip[1,i]) else 
      xylim.clip[2,i] <- min(q$threshold, xylim.clip[2,i])
    xylim.clip
  } 
  if(q$threshold > xylim[1,i]) lines(q$lower, clip(xylim, i, FALSE), ...)
  if(q$threshold < xylim[2,i]) lines(q$upper, clip(xylim, i, TRUE), ...)
  xlim <- xylim[, j]
  xy <- cbind(c(q$threshold, q$threshold), xlim)
  lines(xy[, order(i:j)],  ...)
}
lines.quadtree.leaf <- function(q, xylim, ...) {} # Nothing to do at leaves!

Comme autre exemple, j'ai généré 1 000 000 de points et les ai répartis en groupes de 5 à 9 chacun. Le chronométrage était de 91,7 secondes.

n <- 25000       # Points per cluster
n.centers <- 40  # Number of cluster centers
sd <- 1/2        # Standard deviation of each cluster
set.seed(17)
centers <- matrix(runif(n.centers*2, min=c(-90, 30), max=c(-75, 40)), ncol=2, byrow=TRUE)
xy <- matrix(apply(centers, 1, function(x) rnorm(n*2, mean=x, sd=sd)), ncol=2, byrow=TRUE)
k <- 5
system.time(qt <- quadtree(xy, k))
#
# Set up to map the full extent of the quadtree.
#
xylim <- cbind(x=c(min(xy[,1]), max(xy[,1])), y=c(min(xy[,2]), max(xy[,2])))
plot(xylim, type="n", xlab="x", ylab="y", main="Quadtree")
#
# This is all the code needed for the plot!
#
lines(qt, xylim, col="Gray")
points(qt, pch=".")

entrez la description de l'image ici


Pour illustrer comment interagir avec un SIG , écrivons toutes les cellules quadtree sous forme de fichier de formes polygonales à l'aide de la shapefilesbibliothèque. Le code émule les routines d'écrêtage de lines.quadtree, mais cette fois il doit générer des descriptions vectorielles des cellules. Ils sont sortis sous forme de trames de données à utiliser avec la shapefilesbibliothèque.

cell <- function(q, xylim, ...) {
  if (class(q)=="quadtree") f <- cell.quadtree else f <- cell.quadtree.leaf
  f(q, xylim, ...)
}
cell.quadtree <- function(q, xylim, ...) {
  i <- q$index
  j <- 3 - q$index
  clip <- function(xylim.clip, i, upper) {
    if (upper) xylim.clip[1, i] <- max(q$threshold, xylim.clip[1,i]) else 
      xylim.clip[2,i] <- min(q$threshold, xylim.clip[2,i])
    xylim.clip
  } 
  d <- data.frame(id=NULL, x=NULL, y=NULL)
  if(q$threshold > xylim[1,i]) d <- cell(q$lower, clip(xylim, i, FALSE), ...)
  if(q$threshold < xylim[2,i]) d <- rbind(d, cell(q$upper, clip(xylim, i, TRUE), ...))
  d
}
cell.quadtree.leaf <- function(q, xylim) {
  data.frame(id = q$id, 
             x = c(xylim[1,1], xylim[2,1], xylim[2,1], xylim[1,1], xylim[1,1]),
             y = c(xylim[1,2], xylim[1,2], xylim[2,2], xylim[2,2], xylim[1,2]))
}

Les points eux-mêmes peuvent être lus directement en utilisant read.shpou en important un fichier de données de coordonnées (x, y).

Exemple d'utilisation:

qt <- quadtree(xy, k)
xylim <- cbind(x=c(min(xy[,1]), max(xy[,1])), y=c(min(xy[,2]), max(xy[,2])))
polys <- cell(qt, xylim)
polys.attr <- data.frame(id=unique(polys$id))
library(shapefiles)
polys.shapefile <- convert.to.shapefile(polys, polys.attr, "id", 5)
write.shapefile(polys.shapefile, "f:/temp/quadtree", arcgis=TRUE)

(Utilisez l'étendue souhaitée xylimpour ouvrir une fenêtre dans une sous-région ou étendre le mappage à une région plus grande; ce code correspond par défaut à l'étendue des points.)

Cela suffit à lui seul: une jonction spatiale de ces polygones aux points d'origine identifiera les clusters. Une fois identifiées, les opérations de "résumé" de la base de données généreront des statistiques récapitulatives des points dans chaque cellule.

whuber
la source
Hou la la! Fantastique. Je vais lui donner un coup de feu avec mes données une fois de retour au bureau =)
radek
4
Meilleure réponse @whuber! +1
MerseyViking
1
(1) Vous pouvez lire les fichiers de formes directement avec ( entre autres ) le shapefilespaquet ou bien vous pouvez exporter les coordonnées (x, y) en texte ASCII et les lire avec read.table. (2) Je recommande d'écrire qtsous deux formes: premièrement, en tant que fichier de formes ponctuel pour indiquer xyoù les idchamps sont inclus en tant qu'identificateurs de cluster; deuxièmement, lorsque les segments de ligne tracés par lines.quadtreesont écrits sous forme de fichier de formes polyligne (ou lorsque le traitement analogue écrit les cellules sous forme de fichier de formes polygonales). C'est aussi simple que de modifier lines.quadtree.leafpour produire xylimun rectangle. (Voir les modifications.)
whuber
1
@whubber Merci beaucoup pour une mise à jour. Tout s'est bien passé. Bien mérité +50, bien que maintenant je pense qu'il mérite +500!
radek
1
Je soupçonne que les identifiants calculés n'étaient pas uniques pour une raison quelconque. Apportez ces modifications dans la définition de quad: (1) initialiser id=1; (2) le changement id/2de id*2la lower=ligne; (3) effectuer une modification similaire à id*2+1dans la upper=ligne. (Je vais modifier ma réponse pour refléter cela.) Cela devrait également prendre en charge le calcul de la zone: selon votre SIG, toutes les zones seront positives ou toutes seront négatives. S'ils sont tous négatifs, inversez les listes pour xet ydedans cell.quadtree.leaf.
whuber
6

Voyez si cet algorithme donne suffisamment d'anonymat pour votre échantillon de données:

  1. commencer par une grille régulière
  2. si le polygone a moins que le seuil, fusionnez avec le voisin alternant (E, S, W, N) en spirale dans le sens horaire.
  3. si le polygone a moins que le seuil, passez à 2, sinon passez au polygone suivant

Par exemple, si le seuil minimum est 3:

algorithme

Paulo Scardine
la source
1
Le diable est dans les détails: il semble que cette approche (ou presque n'importe quel algorithme de clustering agglomératif) menace de laisser de petits points "orphelins" dispersés partout, qui ne peuvent alors pas être traités. Je ne dis pas que cette approche est impossible, mais je maintiendrais un scepticisme sain en l'absence d'un algorithme réel et d'un exemple de son application à un ensemble de données ponctuelles réalistes.
whuber
En effet, cette approche pourrait être problématique. Une application de cette méthode à laquelle je pensais utilise des points comme représentations de bâtiments résidentiels. Je pense que cette méthode fonctionnerait bien dans les zones plus densément peuplées. Cependant, il y aurait encore des cas où il y a littéralement un ou deux bâtiments `` au milieu de nulle part '' et cela prendrait beaucoup d'itérations et entraînerait de très grandes surfaces pour enfin atteindre le seuil minimum.
radek
5

De manière similaire à la solution intéressante de Paulo, que diriez-vous d'utiliser un algorithme de subdivision en quatre arbres?

Définissez la profondeur à laquelle vous souhaitez que le quadtree se rende. Vous pouvez également avoir un nombre minimum ou maximum de points par cellule afin que certains nœuds soient plus profonds / plus petits que d'autres.

Subdivisez votre monde, jetez les nœuds vides. Rincer et répéter jusqu'à ce que les critères soient remplis.

MerseyViking
la source
Merci. Quel logiciel recommanderiez-vous pour cela?
radek
1
En principe, c'est une bonne idée. Mais comment surgiraient des nœuds vides si vous n'autorisez jamais moins d'un nombre minimum positif de points par cellule? (Il existe de nombreux types d'arbres quadruples, donc la possibilité de nœuds vides indique que vous en avez un qui n'est pas adapté aux données, ce qui soulève des inquiétudes quant à son utilité pour la tâche prévue.)
whuber
1
Je l'imagine comme ceci: imaginez qu'un nœud a plus que le seuil maximal de points en lui, mais ils sont regroupés vers le haut à gauche du nœud. Le nœud sera subdivisé, mais le nœud enfant en bas à droite sera vide, il peut donc être élagué.
MerseyViking
1
Je vois ce que tu fais (+1). L'astuce consiste à subdiviser en un point déterminé par les coordonnées (telles que leur médiane), garantissant ainsi aucune cellule vide. Sinon, le quadtree est déterminé principalement par l'espace occupé par les points et non par les points eux-mêmes; votre approche devient alors un moyen efficace de concrétiser l'idée générique proposée par @Paulo.
whuber