Menghitung jarak maksimum dalam poligon dalam arah x (arah timur-barat) di PostGIS?

Saya tertarik pada lebar maksimum poligon misalnya danau, ke arah timur-barat. Bounding box hanya akan membantu dalam poligon sederhana tetapi tidak dalam poligon cekung yang kompleks.

Ini kemungkinan membutuhkan beberapa skrip di platform GIS.

Metode yang paling efisien (asimptotik) adalah sapuan garis vertikal: ia membutuhkan penyortiran tepi berdasarkan koordinat y minimumnya dan kemudian memproses tepian dari bawah (minimum y) ke atas (maksimum y), untuk O (e * log) e)) algoritma ketika e edge terlibat.

Prosedurnya, meskipun sederhana, secara mengejutkan sulit untuk dilakukan dengan benar dalam semua kasus. Poligon bisa jahat: mereka dapat memiliki menggantung, sliver, lubang, terputus, memiliki duplikasi simpul, menjalankan simpul di sepanjang garis lurus, dan memiliki batas yang tidak larut antara dua komponen yang berdekatan. Berikut adalah contoh yang menunjukkan banyak karakteristik ini (dan banyak lagi):

Kami akan secara khusus mencari segmen horizontal panjang maksimum yang terletak seluruhnya di dalam penutupan poligon. Sebagai contoh, ini menghilangkan menjuntai antara x = 20 dan x = 40 yang berasal dari lubang antara x = 10 dan x = 25. Kemudian langsung untuk menunjukkan bahwa setidaknya satu dari segmen horisontal dengan panjang maksimum memotong setidaknya satu titik. (Jika ada solusi berpotongan ada simpul, mereka akan berbohong di pedalaman beberapa genjang dibatasi di bagian atas dan bawah dengan solusi yang melakukan berpotongan setidaknya satu titik. Ini memberi kita sarana untuk menemukan semua solusi.)

Oleh karena itu, sapuan garis harus dimulai dengan simpul terendah dan kemudian bergerak ke atas (yaitu, menuju nilai y yang lebih tinggi) untuk berhenti di setiap simpul. Di setiap perhentian, kami menemukan tepi baru yang memancar ke atas dari ketinggian itu; menghilangkan ujung yang berhenti dari bawah pada ketinggian itu (ini adalah salah satu ide kunci: itu menyederhanakan algoritma dan menghilangkan setengah dari pemrosesan potensial); dan hati-hati memproses setiap tepi yang terletak seluruhnya pada ketinggian konstan (tepi horizontal).

Misalnya, pertimbangkan status saat tingkat y = 10 tercapai. Dari kiri ke kanan, kami menemukan tepi berikut:

      x.min x.max y.min y.max
 [1,]    10     0     0    30
 [2,]    10    24    10    20
 [3,]    20    24    10    20
 [4,]    20    40    10    10
 [5,]    40    20    10    10
 [6,]    60     0     5    30
 [7,]    60    60     5    30
 [8,]    60    70     5    20
 [9,]    60    70     5    15
[10,]    90   100    10    40

Dalam tabel ini, (x.min, y.min) adalah koordinat titik akhir bawah tepi dan (x.max, y.max) adalah koordinat titik akhir atasnya. Pada tingkat ini (y = 10), tepi pertama dicegat di dalam interiornya, yang kedua dicegat di bagian bawahnya, dan seterusnya. Beberapa sisi yang berakhir pada level ini, seperti dari (10,0) hingga (10,10), tidak termasuk dalam daftar.

Untuk menentukan di mana titik-titik interior dan yang eksterior, bayangkan mulai dari kiri ekstrim - yang di luar poligon, tentu saja - dan bergerak secara horizontal ke kanan. Setiap kali kami melewati tepi yang tidak horizontal , kami berganti-ganti dari eksterior ke interior dan kembali. (Ini adalah ide kunci lain.) Namun, semua titik dalam setiap tepi horizontal ditentukan berada di dalam poligon, apa pun yang terjadi. (Penutupan poligon selalu menyertakan tepinya.)

Melanjutkan contoh, berikut adalah daftar koordinat x di mana ujung non-horizontal dimulai pada atau melewati garis y = 10:

x.array    6.7 10 20 48 60 63.3 65 90
interior     1  0  1  0  1    0  1  0

(Perhatikan bahwa x = 40 tidak ada dalam daftar ini.) Nilai-nilai interiorarray menandai titik akhir kiri segmen interior: 1 menetapkan interval interior, 0 interval eksterior. Jadi, 1 pertama menunjukkan interval dari x = 6,7 ke x = 10 berada di dalam poligon. Angka 0 berikutnya menunjukkan interval dari x = 10 ke x = 20 berada di luar poligon. Dan begitulah hasilnya: array mengidentifikasi empat interval terpisah sebagai di dalam poligon.

Beberapa interval ini, seperti yang dari x = 60 ke x = 63.3, tidak memotong setiap simpul: pemeriksaan cepat terhadap koordinat x dari semua simpul dengan y = 10 menghilangkan interval tersebut.

Selama pemindaian, kami dapat memantau panjang interval ini, mempertahankan data mengenai interval panjang maksimum yang ditemukan sejauh ini.

Perhatikan beberapa implikasi dari pendekatan ini. Verteks berbentuk "v", ketika ditemui, adalah asal dari dua sisi. Oleh karena itu dua sakelar terjadi ketika melewatinya. Switch itu batal. "V" terbalik mana pun bahkan tidak diproses, karena kedua ujungnya dihilangkan sebelum memulai pemindaian kiri-ke-kanan. Dalam kedua kasus tersebut, simpul semacam itu tidak menghalangi segmen horizontal.

Lebih dari dua sisi dapat berbagi titik: ini diilustrasikan pada (10,0), (60,5), (25, 20), dan - meskipun sulit diketahui - pada (20,10) dan (40) , 10). (Itu karena menjuntai pergi (20,10) -> (40,10) -> (40,0) -> (40, -50) -> (40, 10) -> (20, 10). Perhatikan bagaimana simpul di (40,0) juga di bagian dalam tepi yang lain ... itu tidak menyenangkan.) Algoritma ini menangani situasi itu dengan baik.

Situasi rumit diilustrasikan di bagian paling bawah: koordinat x dari segmen non-horizontal ada

30, 50

Ini menyebabkan semuanya di sebelah kiri x = 30 dianggap eksterior, semuanya antara 30 dan 50 menjadi interior, dan semuanya setelah 50 menjadi eksterior lagi. Vertex di x = 40 bahkan tidak pernah dipertimbangkan dalam algoritma ini.

Beginilah bentuk poligon di akhir pemindaian. Saya menunjukkan semua interval interior yang mengandung simpul dalam abu-abu gelap, setiap interval panjang maksimum berwarna merah, dan warna simpul sesuai dengan koordinat y mereka. Interval maksimum adalah 64 unit.

Satu-satunya perhitungan geometris yang terlibat adalah menghitung di mana ujung-ujungnya memotong garis-garis horizontal: itu adalah interpolasi linier sederhana. Perhitungan juga diperlukan untuk menentukan segmen interior mengandung simpul: ini betweenness penentuan, dengan mudah dihitung dengan beberapa ketidaksetaraan. Kesederhanaan ini membuat algoritma yang kuat dan sesuai baik untuk representasi koordinat integer dan floating point.

Jika koordinat bersifat geografis , maka garis horizontal benar-benar berada di lingkaran garis lintang. Panjangnya tidak sulit untuk dihitung: cukup gandakan panjang Euclidean mereka dengan kosinus lintang mereka (dalam model bola). Karenanya algoritma ini beradaptasi dengan baik terhadap koordinat geografis. (Untuk menangani pembungkusan sekitar +180 sumur meridian, seseorang mungkin perlu terlebih dahulu menemukan kurva dari kutub selatan ke kutub utara yang tidak melewati poligon. Setelah menyatakan kembali semua koordinat x sebagai perpindahan horizontal relatif terhadap itu kurva, algoritma ini akan menemukan segmen horizontal maksimum dengan benar.)

Berikut ini adalah Rkode yang diterapkan untuk melakukan perhitungan dan membuat ilustrasi.

#
# Plotting functions.
#
points.polygon <- function(p, ...) {
  points(p$v, ...)
}
plot.polygon <- function(p, ...) {
  apply(p$e, 1, function(e) lines(matrix(e[c("x.min", "x.max", "y.min", "y.max")], ncol=2), ...))
}
expand <- function(bb, e=1) {
  a <- matrix(c(e, 0, 0, e), ncol=2)
  origin <- apply(bb, 2, mean)
  delta <-  origin %*% a - origin
  t(apply(bb %*% a, 1, function(x) x - delta))
}
#
# Convert polygon to a better data structure.
#
# A polygon class has three attributes:
#   v is an array of vertex coordinates "x" and "y" sorted by increasing y;
#   e is an array of edges from (x.min, y.min) to (x.max, y.max) with y.max >= y.min, sorted by y.min;
#   bb is its rectangular extent (x0,y0), (x1,y1).
#
as.polygon <- function(p) {
  #
  # p is a list of linestrings, each represented as a sequence of 2-vectors 
  # with coordinates in columns "x" and "y". 
  #
  f <- function(p) {
    g <- function(i) {
      v <- p[(i-1):i, ]
      v[order(v[, "y"]), ]
    }
    sapply(2:nrow(p), g)
  }
  vertices <- do.call(rbind, p)
  edges <- t(do.call(cbind, lapply(p, f)))
  colnames(edges) <- c("x.min", "x.max", "y.min", "y.max")
  #
  # Sort by y.min.
  #
  vertices <- vertices[order(vertices[, "y"]), ]
  vertices <- vertices[!duplicated(vertices), ]
  edges <- edges[order(edges[, "y.min"]), ]

  # Maintaining an extent is useful.
  bb <- apply(vertices <- vertices[, c("x","y")], 2, function(z) c(min(z), max(z)))

  # Package the output.
  l <- list(v=vertices, e=edges, bb=bb); class(l) <- "polygon"
  l
}
#
# Compute the maximal horizontal interior segments of a polygon.
#
fetch.x <- function(p) {
  #
  # Update moves the line from the previous level to a new, higher level, changing the
  # state to represent all edges originating or strictly passing through level `y`.
  #
  update <- function(y) {
    if (y > state$level) {
      state$level <<- y
      #
      # Remove edges below the new level from state$current.
      #
      current <- state$current
      current <- current[current[, "y.max"] > y, ]
      #
      # Adjoin edges at this level.
      #
      i <- state$i
      while (i <= nrow(p$e) && p$e[i, "y.min"] <= y) {
        current <- rbind(current, p$e[i, ])
        i <- i+1
      }
      state$i <<- i
      #
      # Sort the current edges by x-coordinate.
      #
      x.coord <- function(e, y) {
        if (e["y.max"] > e["y.min"]) {
          ((y - e["y.min"]) * e["x.max"] + (e["y.max"] - y) * e["x.min"]) / (e["y.max"] - e["y.min"])
        } else {
          min(e["x.min"], e["x.max"])
        }
      }
      if (length(current) > 0) {
        x.array <- apply(current, 1, function(e) x.coord(e, y))
        i.x <- order(x.array)
        current <- current[i.x, ]
        x.array <- x.array[i.x]     
        #
        # Scan and mark each interval as interior or exterior.
        #
        status <- FALSE
        interior <- numeric(length(x.array))
        for (i in 1:length(x.array)) {
          if (current[i, "y.max"] == y) {
            interior[i] <- TRUE
          } else {
            status <- !status
            interior[i] <- status
          }
        }
        #
        # Simplify the data structure by retaining the last value of `interior`
        # within each group of common values of `x.array`.
        #
        interior <- sapply(split(interior, x.array), function(i) rev(i)[1])
        x.array <- sapply(split(x.array, x.array), function(i) i[1])

        print(y)
        print(current)
        print(rbind(x.array, interior))


        markers <- c(1, diff(interior))
        intervals <- x.array[markers != 0]
        #
        # Break into a list structure.
        #
        if (length(intervals) > 1) {
          if (length(intervals) %% 2 == 1) 
            intervals <- intervals[-length(intervals)]
          blocks <- 1:length(intervals) - 1
          blocks <- blocks - (blocks %% 2)
          intervals <- split(intervals, blocks)  
        } else {
          intervals <- list()
        }
      } else {
        intervals <- list()
      }
      #
      # Update the state.
      #
      state$current <<- current
    }
    list(y=y, x=intervals)
  } # Update()

  process <- function(intervals, x, y) {
    # intervals is a list of 2-vectors. Each represents the endpoints of
    # an interior interval of a polygon.
    # x is an array of x-coordinates of vertices.
    #
    # Retains only the intervals containing at least one vertex.
    between <- function(i) {
      1 == max(mapply(function(a,b) a && b, i[1] <= x, x <= i[2]))
    }
    is.good <- lapply(intervals$x, between)
    list(y=y, x=intervals$x[unlist(is.good)])
    #intervals
  }
  #
  # Group the vertices by common y-coordinate.
  #
  vertices.x <- split(p$v[, "x"], p$v[, "y"])
  vertices.y <- lapply(split(p$v[, "y"], p$v[, "y"]), max)
  #
  # The "state" is a collection of segments and an index into edges.
  # It will updated during the vertical line sweep.
  #
  state <- list(level=-Inf, current=c(), i=1, x=c(), interior=c())
  #
  # Sweep vertically from bottom to top, processing the intersection
  # as we go.
  #
  mapply(function(x,y) process(update(y), x, y), vertices.x, vertices.y)
}


scale <- 10
p.raw = list(scale * cbind(x=c(0:10,7,6,0), y=c(3,0,0,-1,-1,-1,0,-0.5,0.75,1,4,1.5,0.5,3)),
             scale *cbind(x=c(1,1,2.4,2,4,4,4,4,2,1), y=c(0,1,2,1,1,0,-0.5,1,1,0)),
             scale *cbind(x=c(6,7,6,6), y=c(.5,2,3,.5)))

#p.raw = list(cbind(x=c(0,2,1,1/2,0), y=c(0,0,2,1,0)))
#p.raw = list(cbind(x=c(0, 35, 100, 65, 0), y=c(0, 50, 100, 50, 0)))

p <- as.polygon(p.raw)

results <- fetch.x(p)
#
# Find the longest.
#
dx <- matrix(unlist(results["x", ]), nrow=2)
length.max <- max(dx[2,] - dx[1,])
#
# Draw pictures.
#
segment.plot <- function(s, length.max, colors,  ...) {
  lapply(s$x, function(x) {
      col <- ifelse (diff(x) >= length.max, colors[1], colors[2])
      lines(x, rep(s$y,2), col=col, ...)
    })
}
gray <- "#f0f0f0"
grayer <- "#d0d0d0"
plot(expand(p$bb, 1.1), type="n", xlab="x", ylab="y", main="After the Scan")
sapply(1:length(p.raw), function(i) polygon(p.raw[[i]], col=c(gray, "White", grayer)[i]))
apply(results, 2, function(s) segment.plot(s, length.max, colors=c("Red", "#b8b8a8"), lwd=4))
plot(p, col="Black", lty=3)
points(p, pch=19, col=round(2 + 2*p$v[, "y"]/scale, 0))
points(p, cex=1.25)

Berikut ini adalah solusi berbasis raster. Ini cepat (saya melakukan semua pekerjaan dari awal hingga selesai dalam 14 menit), tidak memerlukan skrip, hanya membutuhkan beberapa operasi, dan cukup akurat.

Mulailah dengan representasi raster poligon. Yang ini menggunakan kisi 550 baris dan 1200 kolom:

Dalam representasi ini, sel abu-abu (di dalam) memiliki nilai 1 dan semua sel lainnya adalah NoData.

Hitung akumulasi aliran dalam arah barat-timur menggunakan nilai-nilai sel unit untuk grid berat (jumlah "curah hujan"):

Akumulasi rendah gelap, meningkat ke akumulasi tertinggi di kuning cerah.

Zonal maksimum (menggunakan poligon untuk kisi-kisi dan akumulasi aliran untuk nilai-nilai) mengidentifikasi sel di mana aliran terbesar. Untuk menunjukkan ini, saya harus memperbesar ke kanan bawah:

Sel-sel merah menandai ujung akumulasi aliran tertinggi: mereka adalah titik akhir paling kanan dari segmen interior maksimal panjang poligon.

Untuk menemukan segmen-segmen ini, tempatkan semua berat pada sel-sel merah dan jalankan aliran mundur!

Garis merah di dekat bagian bawah menandai dua baris sel: di dalamnya terletak segmen horizontal maksimum-panjang. Gunakan representasi ini apa adanya untuk analisis lebih lanjut atau mengubahnya menjadi bentuk polyline (atau poligon).

Ada beberapa kesalahan diskritisasi yang dibuat dengan representasi raster. Ini dapat dikurangi dengan meningkatkan resolusi, dengan biaya dalam waktu komputasi.

Salah satu aspek yang sangat bagus dari pendekatan ini adalah bahwa biasanya kita menemukan nilai-nilai ekstrem dari hal-hal sebagai bagian dari alur kerja yang lebih besar di mana beberapa tujuan perlu dicapai: penempatan pipa atau lapangan sepak bola, menciptakan penyangga ekologis, dan sebagainya. Prosesnya melibatkan pertukaran. Dengan demikian, garis horizontal paling panjang mungkin bukan bagian dari solusi optimal. Kita mungkin peduli bukan untuk tahu di mana hampir garis terpanjang akan berbohong. Ini sederhana: daripada memilih aliran maksimum zona, pilih semua sel yang dekat dengan zona maksimum. Dalam contoh ini, zonal max sama dengan 744 (jumlah kolom yang direntang oleh segmen interior terpanjang). Sebagai gantinya, mari kita pilih semua sel dalam 5% dari maksimum itu:

Menjalankan aliran dari timur ke barat menghasilkan kumpulan segmen horisontal ini:

Ini adalah peta lokasi di mana batas timur-barat yang tidak terputus adalah 95% atau lebih besar dari batas timur-barat maksimal di mana saja dalam poligon.

Baik. Saya punya ide lain (lebih baik) ( ide-№2 ). Tapi saya kira itu lebih baik direalisasikan sebagai skrip python, bukan sebagai SQL-querry. Sekali lagi di sini adalah kasus umum, bukan hanya EW.

Anda akan membutuhkan kotak pembatas untuk poligon dan azimuth (A) sebagai arah pengukuran Anda. Asumsikan bahwa panjang tepi BBox adalah LA dan LB. Maksimum yang mungkin jarak (MD) dalam poligon adalah: MB = (LA^2 * LB^2)^(1/2), sehingga mencari nilai (V) tidak lebih besar dari MB: V <= MB.

1. Mulai dari titik BBox mana pun membuat garis (LL) dengan panjang MB dan azimuth A.
2. Perpotongan garis LL dengan poligon untuk mendapatkan intersection line (IL)
3. Periksa geometri IL - jika hanya ada dua titik dalam garis IL maka hitung panjangnya. Jika 4 atau lebih - hitung segmen dan ambil bagian yang paling panjang. Null (tidak ada persimpangan sama sekali) - abaikan.
4. Terus buat garis LL lain bergerak dari penghitung titik awal atau searah jarum jam menuju tepi BBox sampai Anda tidak akan berakhir di titik awal (Anda akan melakukan seluruh loop melintasi BBox).
5. Ambil nilai panjang IL terbesar (sebenarnya Anda tidak harus menyimpan semua panjang, Anda dapat menyimpan nilai maksimum 'sejauh ini' sambil mengulang) - itu akan menjadi apa yang Anda cari.

Saya tidak yakin bahwa jawaban Fetzer adalah apa yang ingin Anda lakukan, tetapi st_box2d dapat melakukan pekerjaan itu.

Gagasan SS_Rebelious N ° 1 akan berfungsi dalam banyak kasus tetapi tidak untuk beberapa poligon cekung.

Saya pikir Anda harus membuat lw-line buatan yang titik mengikuti tepian ketika garis verteks dibuat melintasi batas poligon jika ada kemungkinan garis timur-barat.

Untuk ini, Anda bisa mencoba membuat 4 simpul poligon di mana panjang garisnya tinggi, buat poligon P * yang merupakan tumpang tindih sebelumnya dengan Anda poligon asli, dan lihat apakah min (y1) dan maks (y2) tinggalkan beberapa x-line kemungkinan. (di mana y1 adalah himpunan titik antara sudut kiri atas dan sudut kanan atas dan y2 himpunan antara sudut kiri bawah dan kanan bawah dari 4 node poligon Anda). Ini tidak semudah itu saya harap Anda akan menemukan alat psql untuk membantu Anda!

Saya punya ide-№1 ( Edit: untuk kasus umum, tidak hanya arah EW, dan dengan beberapa batasan yang dijelaskan dalam komentar). Saya tidak akan memberikan kode, hanya sebuah konsep. "Arah-x" sebenarnya adalah azimuth, yang dihitung oleh ST_Azimuth. Langkah-langkah yang diusulkan adalah:

1. Ekstrak semua simpul dari poligon sebagai poin.
2. Buat garis di antara setiap pasangan poin.
3. Pilih garis (sebut saja garis-lw) yang berada dalam poligon asli (kita tidak perlu garis yang akan melintasi batas poligon).
4. Temukan jarak dan azimuth untuk setiap lw-line.
5. Pilih jarak terpanjang dari garis-lw di mana azimuth sama dengan dicari-untuk azimuth atau terletak pada beberapa interval (mungkin tidak ada azimuth akan sama persis dengan dicari-untuk azimuth).

Lihatlah pertanyaan saya dan jawaban dari Genius Jahat.

Semoga poligon danau Anda memiliki sejumlah titik, Anda dapat membuat garis pada titik-titik ini dengan azimuth (aspek, arah geo). Pilih panjang garis yang cukup besar (bagian ST_MakePoint), sehingga Anda dapat menghitung garis terpendek antara dua garis paling jauh.

Berikut ini sebuah contoh:

Contoh menunjukkan lebar maksimum poligon. Saya memilih ST_ShortestLine (garis merah) untuk pendekatan ini. ST_MakeLine akan meningkatkan nilai (garis biru) dan titik akhir dari garis (kiri bawah) akan mengenai garis biru poligon. Anda harus menghitung jarak dengan centroid dari garis (bantuan) yang dibuat.

Gagasan untuk poligon tidak beraturan atau cekung untuk pendekatan ini. Mungkin Anda harus memotong poligon dengan raster.