Как я могу создать неправильную сетку, содержащую минимум n точек?

Учитывая большую (~ 1 миллион) выборку неравномерно распределенных точек - возможно ли создать нерегулярную сетку (по размеру, но также может быть неправильной формы, если это возможно?), Которая будет содержать указанное минимальное количество из n точек?

Для меня менее важно, если генерируемые «ячейки» такой сетки содержат ровно n точек или хотя бы n точек.

Мне известны такие решения, как genvecgrid в ArcGIS или Create Grid Layer в QGIS / mmgis, однако все они будут создавать регулярные сетки, которые приведут к выводу с пустыми ячейками (меньшая проблема - я мог бы их просто отбросить) или ячейкам с количеством точек. меньше чем n (большая проблема, так как мне нужно решение для объединения этих ячеек, возможно, с использованием некоторых инструментов отсюда ?).

Я гуглил безрезультатно и открыт как для коммерческих (ArcGIS & extensions), так и для бесплатных (Python, PostGIS, R) решений.

Я вижу, MerseyViking рекомендовал дерево quadtree . Я собирался предложить то же самое, и для того, чтобы объяснить это, вот код и пример. Код написан на, Rно должен легко переноситься, скажем, на Python.

Идея удивительно проста: разделить точки примерно пополам в направлении x, затем рекурсивно разделить две половины вдоль направления y, чередуя направления на каждом уровне, пока больше не требуется разделение.

Поскольку цель состоит в том, чтобы замаскировать фактические местоположения точек, полезно внести некоторую случайность в расщепления . Один быстрый и простой способ сделать это - разделить квантильный набор на небольшую случайную величину от 50%. Таким образом, (а) значения расщепления вряд ли будут совпадать с координатами данных, так что точки будут однозначно попадать в квадранты, созданные разделением, и (б) координаты точек будет невозможно точно восстановить из дерева квадрантов.

Поскольку цель состоит в том, чтобы поддерживать минимальное количество kузлов в каждом листе дерева квадрантов, мы реализуем ограниченную форму дерева квадрантов. Он будет поддерживать (1) точки кластеризации в группы, имеющие от k2 до k1 элементов в каждой и (2) отображение квадрантов.

Этот Rкод создает дерево узлов и конечных листьев, различая их по классам. Маркировка класса ускоряет последующую обработку, такую ​​как построение графиков, показанное ниже. Код использует числовые значения для идентификаторов. Это работает до глубины 52 в дереве (при использовании двойных чисел; если используются целые числа без знака, максимальная глубина составляет 32). Для более глубоких деревьев (что маловероятно в любом приложении, потому что по крайней мере k* 2 ^ 52 балла), идентификаторы должны быть строками.

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

Обратите внимание, что рекурсивная конструкция «разделяй и властвуй» этого алгоритма (и, следовательно, большинства алгоритмов постобработки) означает, что требование времени составляет O (m), а использование ОЗУ - O (n), где m- число ячеек и nколичество очков. mпропорционально nделится на минимальное количество очков на ячейку,k, Это полезно для оценки времени вычислений. Например, если требуется 13 секунд для разделения n = 10 ^ 6 точек на ячейки с 50-99 точками (k = 50), m = 10 ^ 6/50 = 20000. Если вы хотите вместо этого разделить до 5-9 точек на ячейку (k = 5), m в 10 раз больше, поэтому время увеличивается примерно до 130 секунд. (Поскольку процесс разделения набора координат вокруг их середин ускоряется по мере того, как ячейки становятся меньше, фактическое время составило всего 90 секунд.) Чтобы перейти к k = 1 точке на ячейку, потребуется примерно в шесть раз больше времени. еще или девять минут, и мы можем ожидать, что код на самом деле будет немного быстрее.

Прежде чем идти дальше, давайте сгенерируем некоторые интересные нерегулярно расположенные данные и создадим их ограниченное квадродерево (прошедшее время 0,29 секунды):

Вот код для создания этих графиков. Он использует Rполиморфизм: points.quadtreeбудет вызываться всякий раз, когда pointsфункция применяется, например, к quadtreeобъекту. Сила этого очевидна в чрезвычайной простоте функции раскраски точек в соответствии с их кластерным идентификатором:

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

Построение самой сетки немного сложнее, поскольку требует многократного ограничения пороговых значений, используемых для разбиения на четыре дерева, но тот же рекурсивный подход прост и элегантен. Используйте вариант для построения многоугольных представлений квадрантов, если это необходимо.

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!

В качестве другого примера я набрал 1 000 000 очков и разделил их на группы по 5-9 в каждой. Время было 91,7 секунды.

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

В качестве примера того, как взаимодействовать с ГИС , давайте выпишем все ячейки квадродерева в виде шейп-файла многоугольника, используя shapefilesбиблиотеку. Код эмулирует подпрограммы отсечения lines.quadtree, но на этот раз он должен генерировать векторные описания ячеек. Они выводятся как фреймы данных для использования с shapefilesбиблиотекой.

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]))
}

Сами точки могут быть прочитаны непосредственно с помощью read.shpили путем импорта файла данных с координатами (x, y).

Пример использования:

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)

(Здесь можно использовать любой желаемый экстент, чтобы xylimперейти в субрегион или расширить отображение в более крупный регион; этот код по умолчанию соответствует экстенту точек.)

Одного этого достаточно: пространственное соединение этих полигонов с исходными точками идентифицирует кластеры. После идентификации операции «суммирования» базы данных будут генерировать сводную статистику точек в каждой ячейке.

Посмотрите, дает ли этот алгоритм достаточную анонимность для вашего образца данных:

1. начать с обычной сетки
2. если у многоугольника меньше порогового значения, объединить с чередующимся (E, S, W, N) соседним узлом по часовой стрелке.
3. если полигон имеет меньше порога, переходите к 2, иначе переходите к следующему полигону

Например, если минимальный порог равен 3:

Как и в случае с интересным решением Пауло, как насчет использования алгоритма подразделения четырехугольных деревьев?

Установите глубину, на которую вы хотите, чтобы дерево квадратов достигло. Вы также можете иметь минимальное или максимальное количество точек на ячейку, чтобы некоторые узлы были глубже / меньше, чем другие.

Разделите свой мир, отбрасывая пустые узлы. Промыть и повторять до тех пор, пока критерии не будут выполнены.

Другой подход заключается в создании очень тонкой сетки и использовании алгоритма max-p. http://pysal.readthedocs.org/en/v1.7/library/region/maxp.html