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

20

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

Radek
sumber
1
Seberapa "biasa" seharusnya grid itu? Saya ingin tahu apakah Anda dapat melakukan beberapa pengelompokan hierarkis dan kemudian hanya memotong dendrogram untuk memenuhi kebutuhan Anda (walaupun ini mungkin memperluas apa yang akan didefinisikan sebagai konfigurasi spasial biasa). Dokumentasi CrimeStat memiliki beberapa contoh bagus tentang jenis pengelompokan ini .
Andy W
5
Bisakah Anda jelaskan apa yang Anda maksud dengan "grid tidak teratur"? Kedengarannya seperti sebuah oxymoron :-). Lebih penting lagi, apa tujuan latihan ini? Perhatikan juga, bahwa kriteria atau kendala tambahan kemungkinan dibutuhkan: setelah semua, jika Anda menggambar kotak sekitar 1 juta poin, itu dapat dianggap sebagai bagian dari kisi dan itu akan berisi lebih dari n dari mereka. Anda mungkin tidak akan peduli dengan solusi sepele ini: tetapi mengapa tidak, tepatnya?
whuber
@AndyW Terima kasih. Ide bagus & layak ditelusuri. Akan lihat. Ukuran & bentuk 'kisi' sangat penting bagi saya - prioritas (karena privasi data) adalah untuk 'menyembunyikan' n fitur di belakang satu
radek
@whuber Terima kasih juga. Saya setuju - tetapi tidak yakin bagaimana lagi saya bisa menyebutkan partisi seperti itu. Seperti disebutkan di atas - motivasi utama saya adalah privasi data. Memiliki lima lokasi titik (yang tidak dapat saya tunjukkan pada peta akhir) saya ingin mewakili mereka berdasarkan area yang meliputi mereka; dan dapatkan mean / median / dll. nilai untuk itu. Saya setuju bahwa itu mungkin untuk menggambar satu persegi panjang atau cembung lambung mewakili mereka semua - yang akan menjadi perlindungan privasi data pamungkas? ;] Namun - akan lebih berguna untuk mewakilinya dengan bentuk yang melintang, katakanlah 10 fitur. Lalu - saya masih bisa melestarikan pola spasial.
radek
1
IMO memberikan deskripsi Anda, saya akan menggunakan beberapa jenis interpolasi dan menampilkan peta raster (mungkin bandwidth adaptif ukuran minimal N Anda akan cukup untuk memuluskan data). Sejauh CrimeStat, file terbesar yang saya gunakan adalah sekitar 100.000 kasus yang saya percaya (dan pengelompokan pasti akan memakan waktu). Kemungkinan Anda dapat melakukan pra-generalisasi data Anda untuk menyatakannya sebagai lebih sedikit kasus dan masih mendapatkan hasil yang diinginkan untuk apa pun yang Anda inginkan. Ini adalah program yang sangat sederhana, saya sarankan hanya mengambil beberapa menit untuk mencobanya dan lihat sendiri.
Andy W

Jawaban:

26

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

Quadtree

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

masukkan deskripsi gambar di sini


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.

whuber
sumber
Wow! Fantastis. Akan mencobanya dengan data saya sekali lagi di kantor =)
radek
4
Jawaban teratas @whuber! +1
MerseyViking
1
(1) Anda dapat membaca shapefile secara langsung dengan ( antara lain ) shapefilespaket atau Anda dapat mengekspor (x, y) koordinat dalam teks ASCII dan membacanya read.table. (2) Saya merekomendasikan penulisan qtdalam dua bentuk: pertama, sebagai titik pembentuk untuk xytempat idbidang dimasukkan sebagai pengidentifikasi kluster; kedua, di mana segmen garis diplot oleh lines.quadtreedituliskan sebagai polyline shapefile (atau di mana pemrosesan analog menulis sel sebagai polyfon shapefile). Ini sesederhana memodifikasi lines.quadtree.leafoutput xylimsebagai persegi panjang. (Lihat hasil edit.)
whuber
1
@ Whubber Terima kasih banyak untuk pembaruan. Semuanya bekerja dengan lancar. Pantas +50, meskipun sekarang saya pikir itu layak +500!
radek
1
Saya menduga id yang dihitung tidak unik karena beberapa alasan. Buat perubahan ini dalam definisi quad: (1) inisialisasi id=1; (2) perubahan id/2ke id*2dalam lower=line; (3) membuat perubahan serupa ke id*2+1dalam upper=baris. (Saya akan mengedit jawaban saya untuk mencerminkan hal itu.) Itu juga harus memperhatikan perhitungan area: tergantung pada GIS Anda, semua area akan positif atau semuanya akan negatif. Jika semuanya negatif, balikkan daftar untuk xdan ymasuk cell.quadtree.leaf.
whuber
6

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:

algoritma

Paulo Scardine
sumber
1
Iblis ada dalam perincian: tampaknya pendekatan ini (atau hampir semua algoritma pengelompokan aglomeratif) mengancam untuk meninggalkan sedikit "anak yatim" yang tersebar di semua tempat, yang kemudian tidak dapat diproses. Saya tidak mengatakan pendekatan ini tidak mungkin, tetapi saya akan mempertahankan skeptisisme yang sehat dengan tidak adanya algoritma aktual dan contoh penerapannya ke dataset titik realistis.
whuber
Memang pendekatan ini mungkin bermasalah. Salah satu aplikasi dari metode ini yang saya pikirkan tentang penggunaan poin sebagai representasi bangunan tempat tinggal. Saya pikir metode ini akan bekerja dengan baik di daerah yang lebih padat penduduknya. Namun, masih akan ada kasus ketika ada satu atau dua bangunan 'di tengah-tengah dari mana' dan akan membutuhkan banyak iterasi & akan menghasilkan area yang sangat besar untuk akhirnya mencapai ambang minimum.
radek
5

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.

MerseyViking
sumber
Terima kasih. Perangkat lunak apa yang akan Anda rekomendasikan untuk itu?
radek
1
Pada prinsipnya ini adalah ide yang bagus. Tetapi bagaimana node kosong akan muncul jika Anda tidak pernah mengizinkan kurang dari jumlah minimum poin positif per sel? (Ada banyak jenis quadtrees, sehingga kemungkinan node kosong menunjukkan Anda memiliki satu node yang tidak disesuaikan dengan data, yang menimbulkan kekhawatiran tentang kegunaannya untuk tugas yang dimaksud.)
whuber
1
Saya membayangkannya seperti ini: bayangkan sebuah simpul memiliki lebih dari ambang batas maksimum poin di dalamnya, tetapi mereka berkerumun ke kiri atas dari simpul tersebut. Node akan dibagi lagi, tetapi simpul anak kanan bawah akan kosong, sehingga dapat dipangkas.
MerseyViking
1
Saya melihat apa yang Anda lakukan (+1). Caranya adalah dengan membagi pada titik yang ditentukan oleh koordinat (seperti median mereka), sehingga tidak menjamin sel kosong. Kalau tidak, quadtree ditentukan terutama oleh ruang yang ditempati oleh titik dan bukan titik itu sendiri; pendekatan Anda kemudian menjadi cara yang efektif untuk melaksanakan ide generik yang diusulkan oleh @Paulo.
whuber