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

É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).

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):

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=".")

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.

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:

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.

Une autre approche consiste à créer une grille très fine et à utiliser l'algorithme max-p. http://pysal.readthedocs.org/en/v1.7/library/region/maxp.html