Menentukan apakah titik dikelilingi menggunakan pemrosesan raster

Saya berusaha memperbaiki proses vektor / python yang saat ini sangat rumit untuk model bahaya alam. Saat ini kami memiliki skrip panjang yang menghasilkan jarak / garis pancing dari titik tertentu untuk menentukan:

1. jenis poligon yang berpotongan (mis. Hutan, rumput, rawa, dll)
2. jarak ke poligon itu
3. berapa banyak dari garis-garis ini berpotongan poligon, untuk menentukan seberapa 'dikelilingi' itu.

Ada lebih banyak yang terlibat, tetapi itulah intinya. Saya mencoba menemukan cara untuk meningkatkan ini dan saat ini saya bingung pada bagian 3. Idenya adalah untuk menentukan apakah suatu titik benar-benar dikelilingi oleh poligon, dalam katakanlah 200m

Jadi pada gambar terlampir, saya ingin titik A ditandai sebagai berisiko lebih tinggi daripada titik B karena benar-benar dikelilingi oleh poligon saya. Ini diulang sekitar 13 juta poin jadi bukan tugas kecil dan saya lebih suka memiliki permukaan untuk mendapatkan nilai dari, daripada menjalankan skrip kami. Saya pikir harus ada variasi alat hidrologi atau jalur biaya untuk melakukan ini, tetapi saya sepertinya tidak bisa mengatasinya.

Bagaimana saya bisa melakukan ini?

Ada solusi jalur biaya tetapi Anda harus membuat kode sendiri. Berikut ini tampilannya ketika diterapkan pada setiap titik dalam gambar dalam pertanyaan (sedikit kasar untuk mempercepat perhitungan):

Sel-sel hitam adalah bagian dari poligon di sekitarnya. Warna-warnanya, mulai dari oranye terang (pendek) hingga biru (panjang), menunjukkan jarak maksimum (hingga maksimal 50 sel) yang dapat dicapai dengan line-of-sight traversal tanpa memotong sel-sel poligon. (Setiap sel di luar batas gambar ini diperlakukan sebagai bagian dari poligon.)

Mari kita bahas cara efisien untuk melakukan itu menggunakan representasi data raster. Dalam representasi ini, semua sel poligon "sekitarnya" akan memiliki, katakanlah, nilai-nilai bukan nol dan setiap sel yang dapat "dilihat" akan memiliki nilai nol.

Langkah 1: Mempersiapkan struktur data lingkungan

Pertama-tama Anda harus memutuskan apa artinya bagi satu sel untuk memblokir yang lain. Salah satu aturan paling adil yang dapat saya temukan adalah ini: menggunakan koordinat integral untuk baris dan kolom (dan dengan asumsi sel kuadrat), mari kita pertimbangkan sel mana yang dapat memblokir sel (i, j) dari tampilan di titik asal (0,0). Saya menominasikan sel (i ', j') yang paling dekat dengan segmen garis yang menghubungkan (i, j) ke (0,0) di antara semua sel yang koordinatnya paling banyak berbeda dari i dan j paling banyak 1. Karena ini tidak selalu menghasilkan solusi yang unik (misalnya, dengan (i, j) = (1,2) keduanya (0,1) dan (1,1) akan bekerja sama baiknya), beberapa cara untuk menyelesaikan ikatan diperlukan. Akan lebih baik untuk resolusi ikatan ini untuk menghormati simetri lingkungan melingkar di grid: meniadakan koordinat atau beralih koordinat menjaga lingkungan ini. Oleh karena itu kita dapat memutuskan sel mana yang diblokir (i,

Menggambarkan aturan ini adalah kode prototipe berikut ditulis dalam R. Kode ini mengembalikan struktur data yang akan nyaman untuk menentukan "dikelilingi" sel sewenang-wenang dalam kisi.

screen <- function(k=1) {
  #
  # Returns a data structure:
  #   $offset is an array of offsets
  #   $screened is a parallel array of screened offset indexes.
  #   $distance is a parallel array of distances.
  # The first index always corresponds to (0,0).
  #
  screened.by <- function(xy) {
    uv <- abs(xy)
    if (reversed <- uv[2] > uv[1]) {
      uv <- rev(uv)
    }
    i <- which.min(c(uv[1], abs(uv[1]-uv[2]), uv[2]))
    ij <- uv + c(floor((1-i)/3), floor(i/3)-1)
    if (reversed) ij <- rev(ij)
    return(ij * sign(xy))
  }
  #
  # For each lattice point within the circular neighborhood,
  # find the unique lattice point that screens it from the origin.
  #
  xy <- subset(expand.grid(x=(-k:k), y=(-k:k)), 
               subset=(x^2+y^2 <= k^2) & (x != 0 | y != 0))
  g <- t(apply(xy, 1, function(z) c(screened.by(z), z)))
  #
  # Sort by distance from the origin.
  #
  colnames(g) <- c("x", "y", "x.to", "y.to")
  ij <- unique(rbind(g[, 1:2], g[, 3:4]))
  i <- order(abs(ij[,1]), abs(ij[,2])); ij <- ij[i, , drop=FALSE]
  rownames(ij) <- 1:length(i)
  #
  # Invert the "screened by" relation to produce the "screened" relation.
  #
  # (Row, column) offsets.
  ij.df <- data.frame(ij, i=1:length(i))
  #
  # Distances from the origin (in cells).
  distance <- apply(ij, 1, function(u) sqrt(sum(u*u)))
  #
  # "Screens" relation (represented by indexes into ij).
  g <- merge(merge(g, ij.df), ij.df, 
             by.x=c("x.to", "y.to"), by.y=c("x","y"))
  g <- subset(g, select=c(i.x, i.y))
  h <- by(g$i.y, g$i.x, identity)

  return( list(offset=ij, screened=h, distance=distance) )
}

Nilai screen(12)digunakan untuk menghasilkan penggambaran relasi skrining ini: panah menunjuk dari sel ke sel yang segera menyaringnya. Coraknya proporsional berdasarkan jarak ke titik asal, yang berada di tengah lingkungan ini:

Perhitungan ini cepat dan perlu dilakukan hanya sekali untuk lingkungan tertentu. Misalnya, ketika melihat 200 m pada grid dengan sel 5 m, ukuran lingkungan adalah 200/5 = 40 unit.

Langkah 2: Menerapkan perhitungan ke titik yang dipilih

Selebihnya mudah: untuk menentukan apakah sel yang terletak di (x, y) (dalam koordinat baris dan kolom) "dikelilingi" sehubungan dengan struktur data lingkungan ini, lakukan pengujian secara rekursif dimulai dengan offset (i, j) = (0,0) (asal lingkungan). Jika nilai dalam kisi poligon di (x, y) + (i, j) adalah nol, maka visibilitas diblokir di sana. Jika tidak, kita harus mempertimbangkan semua offset yang bisa diblokir pada offset (i, j) (yang ditemukan pada waktu O (1) menggunakan struktur data yang dikembalikan oleh screen). Jika tidak ada yang diblokir, kami telah mencapai batas dan menyimpulkan bahwa (x, y) tidak dikelilingi, jadi kami menghentikan perhitungan (dan tidak repot-repot memeriksa titik yang tersisa di lingkungan).

Kami dapat mengumpulkan informasi yang bahkan lebih berguna dengan melacak jarak garis pandang terjauh yang dicapai selama algoritma. Jika ini kurang dari radius yang diinginkan, sel dikelilingi; selain itu tidak.

Berikut adalah Rprototipe dari algoritma ini. Ini lebih lama daripada yang terlihat, karena Rtidak mendukung struktur stack (sederhana) yang diperlukan untuk mengimplementasikan rekursi, sehingga stack juga harus dikodekan. Algoritme yang sebenarnya dimulai sekitar dua-pertiga dari jalan melalui dan hanya membutuhkan selusin baris atau lebih. (Dan setengah dari mereka hanya menangani situasi di sekitar tepi kisi, memeriksa indeks di luar jangkauan dalam lingkungan. Ini dapat dibuat lebih efisien hanya dengan memperluas kisi poligon dengan kbaris dan kolom di sekelilingnya, menghilangkan semua perlu untuk memeriksa kisaran indeks dengan biaya sedikit lebih banyak RAM untuk memegang grid poligon.)

#
# Test a grid point `ij` for a line-of-sight connection to the perimeter
# of a circular neighborhood.  
#   `xy` is the grid.
#   `counting` determines whether to return max distance or count of stack ops.
#   `perimeter` is the assumed values beyond the extent of `xy`.
#
# Grid values of zero admit light; all others block visibility
# Returns maximum line-of-sight distance found within `nbr`.
#
panvisibility <- function(ij, xy, nbr=screen(), counting=FALSE, perimeter=1) {
  #
  # Implement a stack for the algorithm.
  #
  count <- 0 # Stack count
  stack <- list(ptr=0, s=rep(NA, dim(nbr$offset)[1]))
  push <- function(x) {
    n <- length(x)
    count <<- count+n         # For timing
    stack$s[1:n + stack$ptr] <<- x
    stack$ptr <<- stack$ptr+n
  }
  pop <- function() {
    count <<- count+1         # For timing
    if (stack$ptr <= 0) return(NULL)
    y <- stack$s[stack$ptr]
    #stack$s[stack$ptr] <<- NA # For debugging
    stack$ptr <<- stack$ptr - 1
    return(y)
  }
  #
  # Initialization.
  #
  m <- dim(xy)[1]; n <- dim(xy)[2]
  push(1) # Stack the *indexes* of nbr$offset and nbr$screened.
  dist.max <- -1
  #
  # The algorithm.
  #
  while (!is.null(i <- pop())) {
    cell <- nbr$offset[i, ] + ij
    if (cell[1] <= 0 || cell[1] > m || cell[2] <= 0 || cell[2] > n) {
      value <- perimeter
    } else {  
      value <- xy[cell[1], cell[2]]
    }
    if (value==0) {
      if (nbr$distance[i] > dist.max) dist.max <- nbr$distance[i]
      s <- nbr$screened[[paste(i)]]
      if (is.null(s)) {
        #exited = TRUE
        break
      }
      push(s)
    }
  }
  if (counting) return ( count )
  return(dist.max)
}

Dalam contoh ini, sel-sel poligonal berwarna hitam. Warna memberikan jarak garis pandang maksimum (hingga 50 sel) untuk sel non-poligon, mulai dari oranye terang untuk jarak pendek hingga biru gelap untuk jarak terpanjang. (Sel-selnya satu unit lebar dan tinggi.) Garis-garis yang terlihat jelas dibuat oleh "pulau" poligon kecil di tengah "sungai": masing-masing memblok garis panjang sel-sel lain.

Analisis algoritma

Struktur tumpukan mengimplementasikan pencarian kedalaman-pertama dari grafik visibilitas lingkungan untuk bukti bahwa sel tidak dikelilingi. Di mana sel-sel jauh dari poligon apa pun, pencarian ini akan memerlukan pemeriksaan hanya sel O (k) untuk lingkungan melingkar radius-k. Kasus terburuk terjadi ketika ada sejumlah kecil sel poligon yang tersebar di lingkungan tersebut tetapi meskipun demikian batas dari lingkungan tersebut tidak dapat dijangkau: ini membutuhkan pemeriksaan hampir semua sel di setiap lingkungan, yang merupakan O (k ^ 2) operasi.

Perilaku berikut adalah tipikal dari apa yang akan ditemui. Untuk nilai kecil k, kecuali jika poligon mengisi sebagian besar grid, sebagian besar sel non-poligon akan jelas tidak dikelilingi dan algoritme berskala seperti O (k). Untuk nilai menengah, penskalaan mulai terlihat seperti O (k ^ 2). Ketika k menjadi sangat besar, sebagian besar sel akan dikelilingi dan fakta itu dapat ditentukan dengan baik sebelum seluruh lingkungan diperiksa: upaya komputasi algoritma dengan demikian mencapai batas praktis. Batas ini dicapai ketika jari-jari lingkungan mendekati diameter daerah non-poligon terhubung terbesar di grid.

Sebagai contoh, saya menggunakan countingopsi yang dikodekan ke dalam prototipe screenuntuk mengembalikan jumlah operasi stack yang digunakan dalam setiap panggilan. Ini mengukur upaya komputasi. Grafik berikut memplot jumlah rata-rata stack ops sebagai fungsi dari radius lingkungan. Ini menunjukkan perilaku yang diprediksi.

Kita dapat menggunakan ini untuk memperkirakan perhitungan yang dibutuhkan untuk mengevaluasi 13 juta poin pada grid. Misalkan lingkungan k = 200/5 = 40 digunakan. Kemudian beberapa ratus operasi tumpukan akan dibutuhkan rata-rata (tergantung pada kompleksitas kisi poligon dan di mana 13 juta titik terletak relatif terhadap poligon), yang menyiratkan bahwa dalam bahasa yang dikompilasi secara efisien, paling banyak beberapa ribu operasi numerik sederhana akan diperlukan (menambah, mengalikan, membaca, menulis, mengimbangi, dll). Sebagian besar PC akan dapat mengevaluasi sekitar satu juta poin pada tingkat itu. (ItuRimplementasi jauh lebih lambat daripada itu, karena buruk pada algoritma semacam ini, yang mengapa hanya dapat dianggap sebagai prototipe.) Oleh karena itu, kita mungkin berharap bahwa implementasi yang efisien dalam bahasa yang cukup efisien dan sesuai - C ++ dan Python datang ke pikiran - bisa menyelesaikan evaluasi 13 juta poin dalam satu menit atau kurang, dengan asumsi seluruh kotak poligon berada di RAM.

Ketika kisi terlalu besar untuk masuk ke dalam RAM, prosedur ini dapat diterapkan pada bagian ubin kisi. Mereka hanya perlu tumpang tindih dengan kbaris dan kolom; ambil maxima saat tumpang tindih saat meratapi hasilnya.

Aplikasi lain

" Pengambilan " badan air berhubungan erat dengan "kelengkungan" titik-titiknya. Faktanya, jika kita menggunakan radius lingkungan yang sama dengan atau lebih besar dari diameter badan air, kita akan membuat kisi-kisi (non-directional) yang diambil pada setiap titik dalam badan air tersebut. Dengan menggunakan radius lingkungan yang lebih kecil kita setidaknya akan memperoleh batas bawah untuk pengambilan di semua titik pengambilan tertinggi, yang dalam beberapa aplikasi mungkin cukup baik (dan secara substansial dapat mengurangi upaya komputasi). Varian dari algoritma ini yang membatasi hubungan "disaring oleh" ke arah tertentu akan menjadi salah satu cara untuk menghitung pengambilan secara efisien dalam arah tersebut. Perhatikan bahwa varian tersebut memerlukan modifikasi kode untuk screen; kode untuk panvisibilitytidak berubah sama sekali.

Saya pasti dapat melihat bagaimana orang ingin melakukan ini dengan solusi raster, tetapi bahkan dengan # pengurangan poin, saya akan mengharapkan resolusi yang sangat besar / tinggi dan karenanya sulit untuk memproses kisi atau sekumpulan kisi. Mengingat itu, saya bertanya-tanya apakah mengeksploitasi topologi di gdb mungkin lebih efisien. Anda dapat menemukan semua rongga internal dengan sesuatu seperti:

arcpy.env.workspace = 'myGDB'
arcpy.CreateTopology_management('myGDB', 'myTopology', '')    
arcpy.AddFeatureClassToTopology_management('myTopology', 'myFeatures', '1','1')    
arcpy.AddRuleToTopology_management ('myToplogy', 'Must Not Have Gaps (Area)', 'myFeatures', '', '', '')    
arcpy.ValidateTopology_management('myTopology', 'Full_Extent')
arcpy.ExportTopologyErrors_management('myTopology', 'myGDB', 'topoErrors')
arcpy.FeatureToPolygon_management('topoErrors_line','topoErrorsVoidPolys', '0.1')`

Anda kemudian dapat bekerja dengan topoErrorsVoidPolyspola normal Anda Intersect_analysis()atau apa pun. Anda mungkin perlu dipusingkan dengan mengekstraksi polanya topoErrorsVoidPolys. @whuber memiliki sejumlah posting yang sangat bagus tentang hal-hal semacam ini di tempat lain di gis.stackexchange.com.

Jika Anda benar-benar ingin pergi raster ... Saya akan melakukan sesuatu di sepanjang baris kode semu ini (jangan ngeri hanya karena jelas saya adalah kemunduran AML!: P)

1. rasterize poin ("pts_g") dan polys ("polys_g" (
2. void = regiongroup (con (isnull (polys_g), 1))
3. mungkin perlu melakukan sesuatu untuk memperbaiki rongga untuk menghilangkan poligon eksternal / area semesta terbuka
4. pts_surrounded = con (void, pts_g)

Hanya mengada-ada, jadi mungkin perlu disempurnakan.

Jika Anda menggunakan nilai jarak ambang (di sini Anda berbicara sekitar 200 m), maka solusi terbaik adalah menggunakan analisis vektor:

1) buat buffer 200 m di sekitar setiap titik (berwarna hitam pada ilustrasi)

2) gunakan alat intersect (analisis) antara buffer dan poligon (berwarna biru pada ilustrasi). Ini akan terlihat lebih bagus jika Anda melakukan ini antara batas poligon di sekitarnya dan buffer, tetapi hasil akhirnya adalah sama.

3) gunakan fitur untuk poligon (manajemen) untuk membuat poligon di mana poin Anda benar-benar dikelilingi (berwarna merah pada ilustrasi)

4) pilih lapisan berdasarkan lokasi (manajemen) atau gabungan spasial (analisis) untuk mengidentifikasi titik-titik yang dikelilingi. Penggunaan gabungan spasial memungkinkan Anda memiliki informasi tentang poligon penyematan (area poligon, statistik zona ...) yang dapat berguna untuk pemrosesan lebih lanjut.

Alternatif 2b) Tergantung pada kebutuhan Anda, Anda dapat memilih berdasarkan lokasi poligon di sekitarnya dalam jarak 200 m, Anda kemudian dapat mengidentifikasi beberapa jenis "selungkup" tetapi tidak seketat dalam 2).

Mempertimbangkan "kasus labirin", ini dapat membantu: mengevaluasi berapa lama "melarikan diri" dari lokasi.

• Anda sudah dapat mengecualikan dari poin-poin yang dimasukkan sepenuhnya atau sepenuhnya gratis

• kemudian Anda mengonversi rintangan Anda menjadi raster dan menetapkan nilai ke NoData di mana Anda memiliki poligon, dan ke ukuran sel dalam meteran di mana Anda tidak (ini akan membuat raster biaya Anda).

• ketiga, Anda dapat menghitung jarak biaya menggunakan raster biaya yang baru dibuat

• akhirnya, Anda menggunakan statistik zona sebagai tabel berdasarkan batas penyangga dikonversi menjadi raster (membentuk anulus). Jika Anda dapat melarikan diri ke segala arah, minimumnya harus sekitar 200 (tergantung pada ukuran sel analisis Anda). Tetapi jika Anda berada dalam labirin, maksimum akan lebih besar dari 200. Jadi maksimum statistik zona minus 200 akan menjadi nilai kontinu yang menunjukkan betapa sulitnya untuk "melarikan diri".