Bagaimana cara menghitung sudut di mana dua garis berpotongan di PostGIS?

19

Saya ingin menghitung sudut antara dua garis di mana mereka berpotongan di PostGIS.

Titik awal untuk perhitungan sudut dalam PostGIS tampaknya ST_Azimuth - tetapi itu mengambil poin sebagai input. Pikiran pertama saya adalah mengambil titik akhir dari garis yang berpotongan dan melakukan perhitungan Azimuth pada mereka. Itu tidak cukup baik, karena sebagian besar fitur garis tidak lurus, dan saya tertarik pada sudut di persimpangan. Jadi yang saya hasilkan adalah operasi bersarang yang melewati langkah-langkah berikut:

  1. Identifikasi semua persimpangan antara dua tabel fitur garis.
  2. Buat penyangga yang sangat kecil di sekitar titik persimpangan
  3. Identifikasi titik-titik di mana fitur garis memotong eksterior penyangga (mengambil titik pertama jika ada lebih dari satu - Saya benar-benar hanya tertarik pada apakah sudutnya mendekati 0, 90 atau 180 derajat)
  4. Hitung ST_Azimuth untuk dua poin itu.

SQL lengkapnya agak panjang untuk dikirim di sini, tapi saya sudah mengatasinya di sini jika Anda tertarik. (Ngomong-ngomong, apakah ada cara yang lebih baik daripada membawa semua ladang ke bawah pernyataan WITH?)

Hasilnya tidak terlihat benar, jadi saya jelas melakukan sesuatu yang salah:

contoh keluaran 1 contoh keluaran 2

EDIT Saya redid perhitungan dalam EPSG: 3785 dan hasilnya sedikit berbeda tetapi masih tidak benar:

output di 3785 # 1 output di 3785 # 2

Pertanyaan saya adalah di mana kelemahan dalam proses ini. Apakah saya salah memahami apa yang dilakukan oleh ST_Azimuth? Apakah ada masalah CRS? Sesuatu yang lain sama sekali? Atau mungkin ada cara yang jauh lebih sederhana untuk melakukan ini?

postgis geometry mvexel
sumber
1
Apa CRS asli? Perhitungan sudut harus dibuat dengan proyeksi konformal — bukan dengan lat / long yang tidak terproyeksi (SRID = 4326).
Mike T
Itu EPSG: 4326 koordinat awalnya, saya memasukkan ST_Translate hanya untuk menjadi 100% yakin semua pemrosesan akan dilakukan dalam CRS yang sama. Saya akan mencoba proyeksi konformal, terima kasih.
mvexel
Saya redid perhitungannya adalah EPSG: 3785 dan itu membuat perbedaan - saya akan mengubah pertanyaan untuk menunjukkan hasil baru - tetapi hasilnya masih tidak mencerminkan sudut yang sebenarnya.
mvexel

Jawaban:

12

Saya memiliki pencerahan. Itu agak biasa. Saya meninggalkan satu informasi penting untuk PostGIS untuk menghitung sudut yang tepat.

Yang saya hitung adalah sudut antara hanya dua titik yang memotong eksterior penyangga kecil. Untuk menghitung sudut persimpangan, saya perlu menghitung kedua sudut antara kedua titik di bagian luar buffer dan titik persimpangan dari dua fitur garis dan kurangi.

Saya memperbarui SQL lengkap , tapi inilah yang menonjol:

SELECT
    ...
    abs
    (
        round
        (
            degrees
            (
            ST_Azimuth
            (
                points.point2,
                points.intersection
            )
            -
            ST_Azimuth
            (
                points.point1,
                points.intersection
            )
        )::decimal % 180.0
        ,2
    )
)
AS angle
...
FROM
points
mvexel
sumber
1
Saya sedang berpikir tentang sudut titik buffered pada persimpangan, tetapi saya tidak punya waktu untuk pergi secara detail. Aspek lain adalah unit sudut. Anda perlu mengalikan hasil dalam radian dari ST_Azimuth dengan 180.0 / pi () untuk mendapatkan hasil dalam derajat.
Mike T
Yap terima kasih, saya menggunakan fungsi PostgreSQL derajat () untuk itu.
mvexel
Kujang. (Saya bahkan tidak tahu ada fungsi derajat sampai sekarang.) Akan lebih baik untuk membungkus semua logika ini dalam panggilan fungsi, tetapi saya mengalami kesulitan dalam membuat konsep bagaimana cara kerjanya, yaitu ST_IntersectionAngle(...?
Mike T
Saya benar-benar terkejut bahwa ini bukan fungsi PostGIS. Terima kasih atas tanggapan Anda tentang ini.
mvexel
2

Baru-baru ini saya harus menghitung hal yang sama, tetapi memutuskan pendekatan yang lebih sederhana dan lebih cepat.

Untuk menemukan poin tambahan untuk perhitungan azimuth, saya hanya memeriksa permyriad panjang di belakang persimpangan (atau setelah dalam kasus yang jarang terjadi pada awal baris) menggunakan ST_Line_Locate_Point dan ST_Line_Interpolate_Point :

abs(degrees( 
  ST_Azimuth (
    intersection, 
    ST_Line_Interpolate_Point(
      line1, 
      abs(ST_Line_Locate_Point(line1, intersection) - 0.0001)
    )
  )
  -
  ST_Azimuth (
    intersection, 
    ST_Line_Interpolate_Point(
      line2, 
      abs(ST_Line_Locate_Point(line2, intersection) - 0.0001)
    )
  )
))

Permyriad bersifat arbitrer dan untuk hasil yang lebih konsisten akan lebih baik menggunakan offset absolut. Misalnya, periksa 20m sebelumnya, Anda akan mengubah 0,0001 ke 20/ST_Length(line1)dan 20/ST_Length(line2)masing - masing.

lynxlynxlynx
sumber