Bagaimana cara membuat buffer berorientasi menggunakan arcpy?

Saya ingin membuat buffer berorientasi untuk setiap poligon di shapefile saya menggunakan arcpy. Dengan berorientasi maksud saya bahwa saya memiliki dua sudut a1 dan a2 yang membatasi arah buffer. Ini diwakili dalam grafik di bawah ini:

Ada ide?

Ringkasan

Jawaban ini menempatkan pertanyaan ke dalam konteks yang lebih besar, menjelaskan algoritme yang efisien yang berlaku untuk representasi shapefile fitur (sebagai "vektor" atau "linestrings" poin), menunjukkan beberapa contoh aplikasinya, dan memberikan kode kerja untuk menggunakan atau porting ke dalam lingkungan GIS.

Latar Belakang

Ini adalah contoh pelebaran morfologis. Secara umum penuh, pelebaran "menyebar" titik-titik suatu daerah ke lingkungan mereka; koleksi titik di mana mereka berakhir adalah "pelebaran." Aplikasi dalam SIG sangat banyak: pemodelan penyebaran api, pergerakan peradaban, penyebaran tanaman, dan banyak lagi.

Secara matematis, dan secara umum sangat bagus (tetapi berguna), suatu pelebaran menyebar sekumpulan titik dalam manifold Riemannian (seperti bidang, bola, atau ellipsoid). Penyebaran ditentukan oleh subset dari bundel singgung pada titik-titik ini. Ini berarti bahwa pada masing-masing titik diberikan satu set vektor (arah dan jarak) (saya menyebutnya "lingkungan"); masing-masing vektor ini menggambarkan jalur geodesi yang dimulai pada titik dasarnya. Titik dasar adalah "menyebar" ke ujung semua jalur ini. (Untuk definisi "pelebaran" yang jauh lebih terbatas yang digunakan secara konvensional dalam pemrosesan gambar, lihat artikel Wikipedia . Fungsi penyebaran dikenal sebagai peta eksponensial). dalam geometri diferensial.)

"Buffering" dari fitur adalah salah satu contoh paling sederhana dari pelebaran seperti itu: disk dengan radius konstan (radius buffer) dibuat (setidaknya secara konseptual) di sekitar setiap titik fitur. Persatuan disk ini adalah buffer.

Pertanyaan ini menanyakan perhitungan pelebaran yang sedikit lebih rumit di mana penyebaran hanya boleh terjadi dalam rentang sudut tertentu (yaitu, dalam sektor melingkar). Ini masuk akal hanya untuk fitur yang tidak memiliki permukaan melengkung yang cukup (seperti fitur kecil pada bola atau ellipsoid atau fitur apa pun di pesawat). Ketika kami bekerja di pesawat, juga berarti mengarahkan semua sektor ke arah yang sama. (Namun, jika kami memodelkan penyebaran api oleh angin, kami ingin sektor-sektor diorientasikan dengan angin dan ukurannya mungkin bervariasi dengan kecepatan angin juga: itu adalah salah satu motivasi untuk definisi umum pelebaran yang saya berikan. ) (Pada permukaan melengkung seperti ellipsoid, secara umum tidak mungkin mengarahkan semua sektor ke arah yang "sama".)

Dalam keadaan berikut pelebaran relatif mudah dihitung:

• Fiturnya ada di pesawat (yaitu, kami melebarkan peta fitur dan mudah-mudahan peta itu cukup akurat).

• Pelebaran akan konstan : penyebaran di setiap titik fitur akan terjadi dalam lingkungan yang kongruen dengan orientasi yang sama.

• Lingkungan umum ini cembung. Cembung sangat menyederhanakan dan mempercepat perhitungan.

Pertanyaan ini cocok dalam keadaan khusus seperti itu: ia meminta dilatasi poligon sewenang-wenang oleh sektor melingkar yang asal-usulnya (pusat disk dari mana mereka berasal) terletak di titik dasar. Asalkan sektor-sektor tersebut tidak span lebih dari 180 derajat, mereka akan cembung. (Sektor yang lebih besar selalu dapat dibagi dua menjadi dua sektor cembung; penyatuan dua dilatasi yang lebih kecil akan memberikan hasil yang diinginkan.)

Penerapan

Karena kami melakukan perhitungan Euclidean - melakukan penyebaran di pesawat - kami dapat melebarkan titik hanya dengan menerjemahkan lingkungan pelebaran ke titik itu. (Untuk dapat melakukan ini, lingkungan membutuhkan asalyang akan sesuai dengan titik dasar. Misalnya, asal usul sektor-sektor dalam pertanyaan ini adalah pusat lingkaran dari mana mereka dibentuk. Asal ini terjadi terletak pada batas sektor ini. Dalam operasi penyanggaan GIS standar, lingkungan adalah lingkaran dengan asal pada pusatnya; sekarang asal terletak di bagian dalam lingkaran. Memilih asal bukan masalah besar secara komputasi, karena perubahan asal hanya menggeser seluruh pelebaran, tetapi mungkin merupakan masalah besar dalam hal pemodelan fenomena alam. The sectorfungsi dalam kode di bawah ini menggambarkan bagaimana asal dapat ditentukan.)

Melebarkan segmen garis bisa rumit, tetapi untuk lingkungan cembung kita dapat membuat dilatasi sebagai penyatuan pelebaran dari dua titik akhir bersama dengan jajaran genjang yang dipilih dengan cermat. (Untuk kepentingan ruang, saya tidak akan berhenti sejenak untuk membuktikan pernyataan matematis seperti ini, tetapi mendorong pembaca untuk mencoba bukti mereka sendiri karena ini adalah latihan yang berwawasan luas.) Berikut ini adalah ilustrasi menggunakan tiga sektor (ditunjukkan dengan warna merah muda). Mereka memiliki unit jari-jari dan sudut-sudutnya diberikan dalam judul. Segmen garis itu sendiri memiliki panjang 2, horisontal, dan ditampilkan dalam warna hitam:

Jajaran genjang ditemukan dengan menempatkan titik merah muda yang sejauh mungkin dari segmen dalam arah vertikal saja . Ini memberikan dua poin lebih rendah dan dua poin atas sepanjang garis yang sejajar dengan segmen. Kami hanya harus menggabungkan empat poin ke dalam jajaran genjang (ditunjukkan dengan warna biru). Perhatikan, di sebelah kanan, bagaimana ini masuk akal bahkan ketika sektor itu sendiri hanya segmen garis (dan bukan poligon sejati): di sana, setiap titik pada segmen tersebut telah diterjemahkan ke arah 171 derajat timur utara untuk jarak yang berkisar dari 0 hingga 1. Himpunan titik akhir ini adalah jajaran genjang yang ditunjukkan. Rincian perhitungan ini muncul dalam bufferfungsi yang ditentukan dalam dilate.edgeskode di bawah ini.

Untuk melebarkan polyline , kami membentuk gabungan pelebaran poin dan segmen yang membentuknya. Dua baris terakhir dari dilate.edgesmelakukan loop ini.

Melebarkan poligon membutuhkan termasuk interior poligon bersama dengan pelebaran batasnya. (Pernyataan ini membuat beberapa asumsi tentang lingkungan pelebaran. Salah satunya adalah bahwa semua lingkungan mengandung titik (0,0), yang menjamin poligon termasuk dalam pelebarannya. titik poligon tidak akan melampaui pelebaran titik batas. Ini adalah kasus untuk lingkungan yang konstan.)

Mari kita lihat beberapa contoh bagaimana ini bekerja, pertama dengan nonagon (dipilih untuk mengungkapkan detail) dan kemudian dengan lingkaran (dipilih untuk mencocokkan ilustrasi dalam pertanyaan). Contoh-contoh akan terus menggunakan tiga lingkungan yang sama, tetapi menyusut ke jari-jari 1/3.

Dalam gambar ini interior poligon berwarna abu-abu, titik pelebaran (sektor) berwarna merah muda, dan pelebaran tepi (genjang) berwarna biru.

"Lingkaran" sebenarnya hanya 60-gon, tetapi ia dengan baik mendekati sebuah lingkaran.

Performa

Ketika fitur dasar diwakili oleh titik N dan lingkungan pelebaran oleh titik M, algoritma ini membutuhkan upaya O (N M) . Ini harus ditindaklanjuti dengan menyederhanakan kekacauan simpul dan tepi dalam serikat pekerja, yang dapat membutuhkan upaya O (N M log (N M)): itu adalah sesuatu yang diminta oleh GIS untuk dilakukan; kita tidak harus memprogram itu.

Upaya komputasi dapat ditingkatkan menjadi O (M + N) untuk fitur-fitur dasar cembung (karena Anda dapat mengetahui cara melakukan perjalanan di sekitar batas baru dengan menggabungkan daftar simpul yang menggambarkan batas-batas dari dua bentuk asli) dengan tepat. Ini juga tidak perlu dibersihkan.

Ketika lingkungan pelebaran perlahan - lahan mengubah ukuran dan / atau orientasi saat Anda bergerak di sekitar fitur dasar, pelebaran tepi dapat diperkirakan dari lambung cembung penyatuan pelebaran titik akhir. Jika dua lingkungan pelebaran memiliki titik M1 dan M2, ini dapat ditemukan dengan upaya O (M1 + M2) menggunakan algoritma yang dijelaskan dalam Shamos & Preparata, Computational Geometry . Oleh karena itu, membiarkan K = M1 + M2 + ... + M (N) menjadi jumlah total simpul dalam lingkungan pelebaran N, kita dapat menghitung pelebaran dalam waktu O (K * log (K)).

Mengapa kita ingin mengatasi generalisasi seperti itu jika yang kita inginkan hanyalah buffer sederhana? Untuk fitur besar di bumi, lingkungan pelebaran (seperti disk) yang memiliki ukuran konstan pada kenyataannya mungkin memiliki ukuran yang bervariasi pada peta tempat perhitungan ini dilakukan. Dengan demikian kami memperoleh cara untuk membuat perhitungan yang akurat untuk ellipsoid sambil terus menikmati semua keunggulan geometri Euclidean.

Kode

Contoh-contoh yang dihasilkan dengan Rprototipe ini , yang dapat dengan mudah porting ke bahasa favorit Anda (Python, C ++, dll). Dalam struktur itu paralel dengan analisis yang dilaporkan dalam jawaban ini sehingga tidak perlu penjelasan terpisah. Komentar mengklarifikasi beberapa detail.

(Mungkin menarik untuk dicatat bahwa perhitungan trigonometri hanya digunakan untuk membuat contoh fitur - yang merupakan poligon reguler - dan sektor-sektor. Tidak ada bagian dari perhitungan pelebaran yang memerlukan trigonometri.)

#
# Dilate the vertices of a polygon/polyline by a shape.
#
dilate.points <- function(p, q) {
  # Translate a copy of `q` to each vertex of `p`, resulting in a list of polygons.
  pieces <- apply(p, 1, function(x) list(t(t(q)+x)))
  lapply(pieces, function(z) z[[1]]) # Convert to a list of matrices
}
#
# Dilate the edges of a polygon/polyline `p` by a shape `q`. 
# `p` must have at least two rows.
#
dilate.edges <- function(p, q) {
  i <- matrix(c(0,-1,1,0), 2, 2)       # 90 degree rotation
  e <- apply(rbind(p, p[1,]), 2, diff) # Direction vectors of the edges
  # Dilate a single edge from `x` to `x+v` into a parallelogram
  # bounded by parts of the dilation shape that are at extreme distances
  # from the edge.
  buffer <- function(x, v) {
    y <- q %*% i %*% v # Signed distances orthogonal to the edge
    k <- which.min(y)  # Find smallest distance, then the largest *after* it
    l <- (which.max(c(y[-(1:k)], y[1:k])) + k-1) %% length(y)[1] + 1
    list(rbind(x+q[k,], x+v+q[k,], x+v+q[l,], x+q[l,])) # A parallelogram
  }
  # Apply `buffer` to every edge.
  quads <- apply(cbind(p, e), 1, function(x) buffer(x[1:2], x[3:4]))
  lapply(quads, function(z) z[[1]]) # Convert to a list of matrices
}
#----------------------- (This ends the dilation code.) --------------------------#
#
# Display a polygon and its point and edge dilations.
# NB: In practice we would submit the polygon, its point dilations, and edge 
#     dilations to the GIS to create and simplify their union, producing a single
#     polygon.  We keep the three parts separate here in order to illustrate how
#     that polygon is constructed.
#
display <- function(p, d.points, d.edges, ...) {
  # Create a plotting region covering the extent of the dilated figure.
  x <- c(p[,1], unlist(lapply(c(d.points, d.edges), function(x) x[,1])))
  y <- c(p[,2], unlist(lapply(c(d.points, d.edges), function(x) x[,2])))
  plot(c(min(x),max(x)), c(min(y),max(y)), type="n", asp=1, xlab="x", ylab="y", ...)
  # The polygon itself.
  polygon(p, density=-1, col="#00000040")
  # The dilated points and edges.
  plot.list <- function(l, c) lapply(l, function(p) 
                  polygon(p, density=-1, col=c, border="#00000040"))
  plot.list(d.points, "#ff000020")
  plot.list(d.edges, "#0000ff20")
  invisible(NULL) # Doesn't return anything
}
#
# Create a sector of a circle.
# `n` is the number of vertices to use for approximating its outer arc.
#
sector <- function(radius, arg1, arg2, n=1, origin=c(0,0)) {
  t(cbind(origin, radius*sapply(seq(arg1, arg2, length.out=n), 
                  function(a) c(cos(a), sin(a)))))
}
#
# Create a polygon represented as an array of rows.
#
n.vertices <- 60 # Inscribes an `n.vertices`-gon in the unit circle.
angles <- seq(2*pi, 0, length.out=n.vertices+1)
angles <- angles[-(n.vertices+1)]
polygon.the <- cbind(cos(angles), sin(angles))
if (n.vertices==1) polygon.the <- rbind(polygon.the, polygon.the)
#
# Dilate the polygon in various ways to illustrate.
#
system.time({
  radius <- 1/3
  par(mfrow=c(1,3))
  q <- sector(radius, pi/12, 2*pi/3, n=120)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="-30 to 75 degrees")

  q <- sector(radius, pi/3, 4*pi/3, n=180)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="-150 to 30 degrees")

  q <- sector(radius, -9/20*pi, -9/20*pi)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="171 degrees")
})

Waktu komputasi untuk contoh ini (dari gambar terakhir), dengan N = 60 dan M = 121 (kiri), M = 181 (tengah), dan M = 2 (kanan), adalah seperempat detik. Namun, sebagian besar untuk tampilan. Biasanya, Rkode ini akan menangani sekitar N M = 1,5 juta per detik (mengambil hanya 0,002 detik atau lebih untuk melakukan semua contoh perhitungan yang ditunjukkan). Meskipun demikian, penampilan produk M N menyiratkan pelebaran banyak figur atau figur rumit melalui lingkungan terperinci mungkin membutuhkan waktu yang lama, jadi waspadalah! Benchmark timing untuk masalah yang lebih kecil sebelum menangani yang besar. Dalam keadaan seperti itu, orang mungkin mencari solusi berbasis raster (yang jauh lebih mudah diimplementasikan, yang pada dasarnya hanya memerlukan satu perhitungan lingkungan).

Ini cukup luas, tetapi Anda dapat:

1. buffer poligon asli
2. menemukan titik asal sinar "berorientasi" yang akan dibuat pada batas poligon (semacam titik singgung?)
3. buat / perpanjang garis keluar dari titik itu ke jarak di luar buffer menggunakan sudut yang dibahas dalam komentar untuk pertanyaan.
4. memotong garis itu dengan buffer dan poligon asli. Ini mungkin bisa dilakukan bersamaan dengan 3) dengan argumen yang tepat untuk Diperpanjang.
5. ekstrak poligon "berorientasi penyangga" baru dari kumpulan poligon yang dihasilkan