Comment calculer l'angle auquel deux lignes se croisent dans PostGIS?

Je veux calculer l'angle entre deux lignes où elles se croisent dans PostGIS.

Le point de départ des calculs d'angle dans PostGIS semble être ST_Azimuth - mais cela prend des points en entrée. Ma première pensée a été de prendre les extrémités des lignes qui se croisent et d'effectuer un calcul d'azimut sur celles-ci. Ce n'est pas suffisant, car la plupart des entités linéaires ne sont pas droites et je m'intéresse à l'angle à l'intersection. Donc, ce que j'ai trouvé est une opération imbriquée qui passe par les étapes suivantes:

1. Identifiez toutes les intersections entre les deux tables d'entités linéaires.
2. Créer un très petit tampon autour du point d'intersection
3. Identifiez les points où les entités linéaires coupent l'extérieur du tampon (en prenant le premier point s'il y en a plusieurs - je ne suis vraiment intéressé que si l'angle est proche de 0, 90 ou 180 degrés)
4. Calculez ST_Azimuth pour ces deux points.

Le SQL complet est un peu long à publier ici, mais je l'ai ici si vous êtes intéressé. (Soit dit en passant, y a-t-il une meilleure façon que de reporter tous les champs en descendant les instructions WITH?)

Les résultats ne semblent pas corrects, donc je fais clairement quelque chose de mal:

EDIT J'ai refait les calculs en EPSG: 3785 et les résultats sont un peu différents mais toujours pas corrects:

Ma question est de savoir où sont les défauts dans ce processus. Suis-je en train de mal comprendre ce que fait ST_Azimuth? Y a-t-il un problème CRS? Quelque chose d'autre? Ou peut-être existe-t-il une manière beaucoup, beaucoup plus simple de le faire?

J'ai eu l'épiphanie. C'est plutôt banal. Je laissais de côté une information essentielle à PostGIS pour calculer l'angle droit.

Ce que je calculais, c'était l'angle entre seulement les deux points coupant l'extérieur du petit tampon. Pour calculer l'angle de l'intersection, je dois calculer les deux angles entre les deux points à l'extérieur du tampon et le point d'intersection des deux entités linéaires et les soustraire.

J'ai mis à jour le SQL complet , mais voici le bit saillant:

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

J'ai récemment dû calculer la même chose, mais j'ai décidé d'une approche plus simple et probablement plus rapide.

Pour trouver les points supplémentaires pour le calcul de l'azimut, je vérifie simplement une permyriade de la longueur derrière l'intersection (ou après dans le cas rare où cela se produit au tout début de la ligne) en utilisant ST_Line_Locate_Point et 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)
    )
  )
))

La permyriade était arbitraire et pour des résultats plus cohérents, il serait préférable d'utiliser un décalage absolu. Pour vérifier par exemple 20 m à l'avance, vous devez changer 0,0001 respectivement 20/ST_Length(line1)et 20/ST_Length(line2).