Jak obliczyć kąt, pod którym przecinają się dwie linie w PostGIS?

Chcę obliczyć kąt między dwiema liniami, w których przecinają się one w PostGIS.

Punktem wyjścia dla obliczeń kąta w PostGIS wydaje się być ST_Azimuth - ale to bierze punkty jako dane wejściowe. Moją pierwszą myślą było wybranie punktów końcowych przecinających się linii i wykonanie na nich obliczeń azymutu. To nie wystarcza, ponieważ większość cech linii nie jest prosta, a mnie interesuje kąt na przecięciu. Tak więc wymyśliłem zagnieżdżoną operację, która przechodzi przez następujące kroki:

1. Zidentyfikuj wszystkie przecięcia między dwiema tabelami cech linii.
2. Utwórz bardzo mały bufor wokół punktu przecięcia
3. Zidentyfikuj punkty, w których cechy linii przecinają zewnętrzną część bufora (biorąc pierwszy punkt, jeśli jest więcej niż jeden - naprawdę interesuje mnie tylko, czy kąt jest bliski 0, 90 lub 180 stopni)
4. Oblicz ST_Azimuth dla tych dwóch punktów.

Pełny SQL jest dość długi, aby opublikować tutaj, ale przedstawiłem go tutaj, jeśli jesteś zainteresowany. (Nawiasem mówiąc, czy istnieje lepszy sposób niż przeniesienie wszystkich pól w dół instrukcji WITH?)

Wyniki nie wyglądają dobrze, więc najwyraźniej robię coś złego:

EDYCJA Zmieniłem obliczenia w EPSG: 3785, a wyniki są trochę inne, ale nadal nie tak:

Moje pytanie brzmi, gdzie są wady tego procesu. Czy nie rozumiem, co robi ST_Azimuth? Czy występuje problem z CRS? Coś jeszcze w ogóle? A może jest na to znacznie prostszy sposób?

Miałem objawienie. To jest raczej przyziemne. Pominąłem jedną istotną informację dla PostGIS, aby obliczyć kąt prosty.

To, co obliczałem, to kąt między tylko dwoma punktami przecinającymi zewnętrzną część bufora. Aby obliczyć kąt przecięcia, muszę obliczyć oba kąty między obydwoma punktami na zewnętrznej stronie zderzaka i punktem przecięcia dwóch elementów linii i odjąć je.

Zaktualizowałem pełny SQL , ale oto najistotniejszy bit:

SELECT
    ...
    abs
    (
        round
        (
            degrees
            (
            ST_Azimuth
            (
                points.point2,
                points.intersection
            )
            -
            ST_Azimuth
            (
                points.point1,
                points.intersection
            )
        )::decimal % 180.0
        ,2
    )
)
AS angle
...
FROM
points 

Niedawno musiałem obliczyć to samo, ale zdecydowałem się na prostsze i prawdopodobnie szybsze podejście.

Aby znaleźć dodatkowe punkty do obliczenia azymutu, po prostu sprawdzam permyriad długości za skrzyżowaniem (lub po rzadkim przypadku, który zdarza się na samym początku linii) za pomocą ST_Line_Locate_Point i 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)
    )
  )
))

Permyriada była arbitralna i dla bardziej spójnych wyników lepiej byłoby zastosować przesunięcie bezwzględne. Aby na przykład wcześniej sprawdzić 20 m, zmień odpowiednio 0,0001 na 20/ST_Length(line1)i 20/ST_Length(line2).