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

20

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

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

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

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

Радек
источник
1
Насколько «регулярной» должна быть сетка? Интересно, можете ли вы сделать некоторую иерархическую кластеризацию, а затем просто вырезать дендрограмму, чтобы удовлетворить ваши потребности (хотя это, вероятно, растягивает то, что можно определить как обычную пространственную конфигурацию). В документации CrimeStat есть несколько хороших примеров кластеризации такого типа .
Энди Ш
5
Не могли бы вы объяснить, что именно вы подразумеваете под "неправильной сеткой"? Это звучит как оксюморон :-). Более конкретно, какова будет цель этого упражнения? Также обратите внимание, что дополнительные критерии или ограничения, вероятно, необходимы: в конце концов, если вы нарисовали квадрат вокруг всех 1 миллиона точек, его можно рассматривать как часть сетки, и он будет содержать более n из них. Вы, вероятно, не заботитесь об этом тривиальном решении, хотя: но почему бы и нет?
whuber
@ AndyW Спасибо. Хорошая идея и стоит изучить. Посмотрим. Размер и форма «сетки» имеют для меня второстепенное значение - приоритет (из-за конфиденциальности данных) состоит в том, чтобы «скрыть» n функций за одним
radek
@whuber Спасибо тоже. Я согласен - но не был уверен, как еще я мог бы назвать такое разбиение. Как уже упоминалось выше, моя главная мотивация - конфиденциальность данных. Имея пять точечных местоположений (которые я не могу показать на итоговой карте), я хотел бы представить их областью, покрывающей их; и получить среднее значение / медиана / и т. д. значение для этого. Я согласен, что было бы возможно нарисовать один прямоугольник или выпуклую оболочку, представляющую их все - это была бы максимальная защита конфиденциальности данных, я полагаю? ;] Однако - было бы более полезно представить его в виде фигур, скажем, 10 объектов. Тогда - я все еще могу сохранить пространственную картину.
Радек
1
IMO, учитывая ваше описание, я использовал бы некоторый тип интерполяции и отобразил растровую карту (возможно, адаптивной полосы пропускания размер вашего минимального N был бы достаточен для сглаживания данных). Что касается CrimeStat, то я считаю, что самые большие файлы, которые я использовал, были около 100 000 случаев (и кластеризация, безусловно, потребовала бы времени). Вполне вероятно, что вы могли бы сделать некоторое предварительное обобщение ваших данных, чтобы представить их как меньшее количество случаев и при этом получить желаемые результаты для всего, что вы хотите. Это действительно простая программа, я бы посоветовал потратить несколько минут, чтобы попробовать и убедиться в этом.
Энди Ш

Ответы:

26

Я вижу, 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перейти в субрегион или расширить отображение в более крупный регион; этот код по умолчанию соответствует экстенту точек.)

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

Whuber
источник
Вот это да! Фантастический.
Сделаю
4
Топ ответ @whuber! +1
MerseyViking
1
(1) Вы можете читать шейп-файлы напрямую из пакета ( среди прочего ) shapefilesили экспортировать (x, y) координаты в текст ASCII и читать их с помощью read.table. (2) Я рекомендую выписать qtв двух формах: во- первых, в качестве точки шейп для xyгде idполя включены в качестве идентификаторов кластера; во-вторых, когда сегменты линии, нанесенные на график lines.quadtree, записываются как шейп-файл полилинии (или когда аналогичная обработка записывает ячейки как шейп-файл многоугольника). Это так же просто, как изменить lines.quadtree.leafвывод xylimв виде прямоугольника. (См.
Правки
1
@ whubber Большое спасибо за обновление. Все работало без сбоев. Хорошо заслужил +50, хотя теперь я думаю, что заслуживает +500!
Радек
1
Я подозреваю, что рассчитанные идентификаторы не были уникальными по некоторым причинам. Внесите эти изменения в определение quad: (1) инициализируйте id=1; (2) изменение , id/2чтобы id*2в lower=линии; (3) внести аналогичные изменения id*2+1в upper=строке. (Я отредактирую свой ответ, чтобы отразить это.) Это также должно учитывать расчет площади: в зависимости от вашей ГИС все области будут положительными или все будут отрицательными. Если они все отрицательные, переверните списки для xи yв cell.quadtree.leaf.
whuber
6

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

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

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

алгоритм

Пауло Скардин
источник
1
Дьявол кроется в деталях: кажется, что этот подход (или почти любой алгоритм агломерационной кластеризации) грозит оставить маленькие «лишние» точки разбросанными повсюду, которые затем не могут быть обработаны. Я не говорю, что такой подход невозможен, но я бы поддерживал здоровый скептицизм в отсутствие реального алгоритма и примера его применения к реалистичному набору точечных данных.
whuber
Действительно, такой подход может быть проблематичным. Одно из применений этого метода, о котором я думал, использует точки как изображения жилых зданий. Я думаю, что этот метод будет хорошо работать в более густонаселенных районах. Однако бывают случаи, когда буквально одно или два здания находятся «посреди ниоткуда», и потребуется много итераций, что приведет к действительно большим площадям, чтобы в итоге достичь минимального порога.
Радек
5

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

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

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

MerseyViking
источник
Благодарю. Какое программное обеспечение вы бы порекомендовали для этого?
Радек
1
В принципе это хорошая идея. Но как возникнут пустые узлы, если вы никогда не допустите меньше положительного минимального количества точек на ячейку? (Существует много видов квадродерев, поэтому наличие пустых узлов означает, что вы имеете в виду тот, который не адаптирован к данным, что вызывает опасения относительно его полезности для намеченной задачи.)
whuber
1
Я представляю это так: представьте, что узел имеет больше, чем максимальный порог точек в нем, но они сгруппированы по направлению к верхнему левому углу узла. Узел будет подразделен, но нижний правый дочерний узел будет пустым, поэтому его можно будет удалить.
MerseyViking
1
Я вижу, что ты делаешь (+1). Хитрость заключается в том, чтобы подразделить в точке, определяемой координатами (например, их медианой), что гарантирует отсутствие пустых ячеек. В противном случае квадродерево определяется в первую очередь пространством, занимаемым точками, а не самими точками; Ваш подход становится эффективным способом реализации общей идеи, предложенной @Paulo.
whuber