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

19

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:

exemple de sortie 1 exemple de sortie 2

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

sortie en 3785 # 1 sortie en 3785 # 2

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?

mvexel
la source
1
Quel était le CRS d'origine? Les calculs d'angle doivent être effectués avec une projection conforme - et non avec un lat / long non projeté (SRID = 4326).
Mike T
C'était des coordonnées EPSG: 4326 à l'origine, j'ai inclus le ST_Translate juste pour être sûr à 100% que tout le traitement se ferait dans le même CRS. Je vais essayer une projection conforme, merci.
mvexel
J'ai refait les calculs en EPSG: 3785 et cela fait une différence - je vais modifier la question pour afficher les nouveaux résultats - mais le résultat ne reflète toujours pas l'angle réel.
mvexel

Réponses:

12

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 
mvexel
la source
1
Je pensais à l'angle du point tamponné par rapport à l'instersection, mais je n'ai pas le temps d'entrer dans les détails. Un autre aspect est les unités angulaires. Vous devez multiplier le résultat en radians de ST_Azimuth par 180,0 / pi () pour obtenir des résultats en degrés.
Mike T
Oui, merci, j'utilise la fonction degrés () de PostgreSQL pour cela.
mvexel
Couperet. (Je ne savais même pas qu'il y avait une fonction degrés jusqu'à présent.) Ce serait bien de résumer toute cette logique dans un appel de fonction, mais j'ai du mal à conceptualiser comment cela fonctionnerait, c'est ST_IntersectionAngle(...-à- dire ?
Mike T
J'ai été en fait surpris que ce ne soit pas une fonction PostGIS. Merci pour vos commentaires à ce sujet.
mvexel
2

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).

lynxlynxlynx
la source