Bagaimana saya bisa menghasilkan grid tidak teratur yang mengandung minimum n poin?

Diberikan sampel besar (~ 1 juta) dari poin yang tidak terdistribusi secara merata - apakah mungkin untuk menghasilkan grid tidak teratur (dalam ukuran, tetapi juga bisa berbentuk tidak teratur jika itu mungkin?) Yang akan berisi jumlah minimum n poin yang ditentukan ?

Itu kurang penting bagi saya jika 'sel' yang dihasilkan dari kisi-kisi tersebut mengandung tepat n jumlah poin atau setidaknya n poin.

Saya mengetahui solusi seperti genvecgrid di ArcGIS atau Create Grid Layer di QGIS / mmgis namun mereka semua akan membuat grid biasa yang akan menghasilkan output dengan sel kosong (masalah yang lebih kecil - saya bisa dengan mudah membuangnya) atau sel dengan jumlah poin kurang dari n (masalah lebih besar karena saya membutuhkan solusi untuk menggabungkan sel-sel tersebut, mungkin menggunakan beberapa alat dari sini ?).

Saya sudah googling di sekitar tetapi tidak berhasil dan saya terbuka untuk solusi komersial (ArcGIS & ekstensi) atau gratis (Python, PostGIS, R).

Saya melihat MerseyViking merekomendasikan quadtree . Saya akan menyarankan hal yang sama dan untuk menjelaskannya, inilah kode dan contohnya. Kode ditulis Rtetapi harus dengan mudah port ke, katakanlah, Python.

Idenya sangat sederhana: pisahkan titik kira-kira setengah dalam arah x, kemudian secara rekursif membagi dua bagian sepanjang arah y, bergantian arah di setiap level, sampai tidak ada lagi pemisahan yang diinginkan.

Karena tujuannya adalah untuk menyamarkan lokasi titik yang sebenarnya, ada baiknya untuk memasukkan beberapa keacakan ke dalam pemisahan . Salah satu cara sederhana cepat untuk melakukan ini adalah dengan membagi pada set kuantil jumlah acak kecil menjauh dari 50%. Dengan cara ini (a) nilai-nilai pemisahan sangat tidak mungkin bertepatan dengan koordinat data, sehingga poin akan jatuh secara unik ke dalam kuadran yang dibuat oleh partisi, dan (b) koordinat titik tidak mungkin untuk merekonstruksi secara tepat dari quadtree.

Karena tujuannya adalah untuk mempertahankan jumlah minimum k node dalam setiap daun quadtree, kami menerapkan bentuk quadtree terbatas. Ini akan mendukung (1) pengelompokan poin ke dalam kelompok yang memiliki antara kdan 2 * k-1 elemen masing-masing dan (2) memetakan kuadran.

Ini RKode menciptakan pohon node dan terminal daun, membedakannya berdasarkan kelas. Pelabelan kelas mempercepat pasca pemrosesan seperti memplot, ditunjukkan di bawah ini. Kode menggunakan nilai numerik untuk id. Ini berfungsi hingga kedalaman 52 di pohon (menggunakan ganda; jika bilangan bulat panjang yang tidak ditandatangani digunakan, kedalaman maksimum adalah 32). Untuk pohon yang lebih dalam (yang sangat tidak mungkin dalam aplikasi apa pun, karena setidaknya k* 2 ^ 52 poin akan terlibat), id harus berupa string.

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

Perhatikan bahwa desain pembagian dan penaklukan rekursif dari algoritma ini (dan, akibatnya, sebagian besar algoritma pasca-pemrosesan) berarti bahwa persyaratan waktu adalah O (m) dan penggunaan RAM adalah O (n) di mana mjumlah sel dan njumlah titik. msebanding dengan ndibagi dengan poin minimum per sel,k. Ini berguna untuk memperkirakan waktu perhitungan. Misalnya, jika dibutuhkan 13 detik untuk mempartisi n = 10 ^ 6 poin ke dalam sel 50-99 poin (k = 50), m = 10 ^ 6/50 = 20000. Jika Anda ingin mempartisi ke 5-9 poin per sel (k = 5), m adalah 10 kali lebih besar, sehingga waktunya naik sekitar 130 detik. (Karena proses pemisahan seperangkat koordinat di sekitar middle mereka semakin cepat karena sel semakin kecil, waktu yang sebenarnya hanya 90 detik.) Untuk mencapai k = 1 titik per sel, akan memakan waktu sekitar enam kali lebih lama masih, atau sembilan menit, dan kita dapat mengharapkan kode sebenarnya menjadi sedikit lebih cepat dari itu.

Sebelum melangkah lebih jauh, mari kita buat beberapa data spasi tidak teratur yang menarik dan buat quadtree terbatas mereka (waktu berlalu 0,29 detik):

Berikut kode untuk membuat plot ini. Ini mengeksploitasi Rpolimorfisme: points.quadtreeakan dipanggil setiap kali pointsfungsi diterapkan ke quadtreeobjek, misalnya. Kekuatan ini terbukti dalam kesederhanaan ekstrim dari fungsi untuk mewarnai titik-titik sesuai dengan pengidentifikasi kluster mereka:

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

Memplot kisi-kisi itu sendiri sedikit lebih sulit karena memerlukan kliping berulang dari ambang yang digunakan untuk partisi quadtree, tetapi pendekatan rekursif yang sama sederhana dan elegan. Gunakan varian untuk membangun representasi poligonal kuadran jika diinginkan.

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!

Sebagai contoh lain, saya menghasilkan 1.000.000 poin dan mempartisinya menjadi 5-9 kelompok. Waktu adalah 91,7 detik.

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

Sebagai contoh bagaimana berinteraksi dengan GIS , mari kita tuliskan semua sel quadtree sebagai bentuk poligon menggunakan shapefilesperpustakaan. Kode ini mengemulasi rutinitas kliping lines.quadtree, tetapi kali ini harus menghasilkan deskripsi vektor sel. Ini adalah output sebagai frame data untuk digunakan dengan shapefilesperpustakaan.

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

Poin itu sendiri dapat dibaca secara langsung menggunakan read.shpatau dengan mengimpor file data koordinat (x, y).

Contoh penggunaan:

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)

(Gunakan batas yang diinginkan untuk di xylimsini untuk membuka jendela ke subkawasan atau untuk memperluas pemetaan ke wilayah yang lebih besar; kode ini default ke tingkat poin.)

Ini saja sudah cukup: gabungan spasial dari poligon-poligon ini ke titik-titik awal akan mengidentifikasi kelompok-kelompok tersebut. Setelah diidentifikasi, operasi "ringkasan" basis data akan menghasilkan ringkasan statistik dari titik-titik dalam setiap sel.

Lihat apakah algoritma ini memberikan anonimitas yang cukup untuk sampel data Anda:

1. mulai dengan kisi-kisi biasa
2. jika poligon memiliki kurang dari ambang batas, gabungkan dengan tetangga yang bergantian (E, S, W, N) berputar searah jarum jam.
3. jika poligon memiliki kurang dari ambang batas, lanjutkan ke 2, lalu pergi ke poligon berikutnya

Misalnya, jika ambang minimum adalah 3:

Demikian pula dengan solusi menarik Paulo, bagaimana dengan menggunakan algoritma subdivisi quad tree?

Atur kedalaman yang Anda inginkan untuk quadtree. Anda juga bisa memiliki jumlah poin minimum atau maksimum per sel sehingga beberapa node akan lebih dalam / lebih kecil dari yang lain.

Bagi dunia Anda, buang node kosong. Bilas dan ulangi sampai kriteria terpenuhi.

Pendekatan lain adalah membuat kisi yang sangat halus, dan menggunakan algoritma max-p. http://pysal.readthedocs.org/en/v1.7/library/region/maxp.html