¿Cómo calcular el ángulo en el que dos líneas se cruzan en PostGIS?

19

Quiero calcular el ángulo entre dos líneas donde se cruzan en PostGIS.

El punto de partida para los cálculos de ángulo en PostGIS parece ser ST_Azimuth , pero eso toma puntos como entrada. Mi primer pensamiento fue tomar los puntos finales de las líneas que se cruzan y realizar un cálculo de acimut en ellos. Eso no es lo suficientemente bueno, porque la mayoría de las entidades de línea no son rectas, y estoy interesado en el ángulo en la intersección. Entonces, se me ocurrió una operación anidada que sigue los siguientes pasos:

  1. Identifique todas las intersecciones entre las dos tablas de entidades de línea.
  2. Cree un búfer muy pequeño alrededor del punto de intersección.
  3. Identifique los puntos donde las entidades de línea se cruzan con el exterior del búfer (tomando el primer punto si hay más de uno; en realidad solo me interesa saber si el ángulo está cerca de 0, 90 o 180 grados)
  4. Calcule ST_Azimuth para esos dos puntos.

El SQL completo es un poco largo para publicar aquí, pero lo mencioné aquí si está interesado. (Por cierto, ¿hay una mejor manera que transferir todos los campos que van por las declaraciones WITH?)

Los resultados no se ven bien, así que claramente estoy haciendo algo mal:

ejemplo de salida 1 ejemplo de salida 2

EDITAR Rehice los cálculos en EPSG: 3785 y los resultados son un poco diferentes, pero aún no son correctos:

salida en 3785 # 1 salida en 3785 # 2

Mi pregunta es dónde están los defectos en este proceso. ¿Estoy malinterpretando lo que hace ST_Azimuth? ¿Hay un problema de CRS? ¿Algo más por completo? ¿O tal vez hay una manera mucho, mucho más simple de hacer esto?

mvexel
fuente
1
¿Cuál fue el CRS original? Los cálculos de ángulo deben hacerse con una proyección conforme, no con lat / long no proyectados (SRID = 4326).
Mike T
Era EPSG: 4326 coordenadas originalmente, incluí el ST_Translate solo para estar 100% seguro de que todo el procesamiento se realizaría en el mismo CRS. Intentaré una proyección conforme, gracias.
mvexel
Rehice los cálculos es EPSG: 3785 y hace la diferencia, enmendaré la pregunta para mostrar los nuevos resultados, pero el resultado aún no refleja el ángulo real.
mvexel

Respuestas:

12

Tuve la epifanía. Es bastante mundano. Estaba dejando de lado una información esencial para que PostGIS calcule el ángulo correcto.

Lo que estaba calculando era el ángulo entre solo los dos puntos que intersectan el pequeño exterior del búfer. Para calcular el ángulo de la intersección, necesito calcular ambos ángulos entre ambos puntos en el exterior del búfer y el punto de intersección de las dos entidades de línea y restarlas.

Actualicé el SQL completo , pero aquí está el bit sobresaliente:

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
fuente
1
Estaba pensando en el ángulo del punto amortiguado con la sección, pero no tengo tiempo para entrar en detalles. Otro aspecto son las unidades angulares. Debe multiplicar el resultado en radianes de ST_Azimuth por 180.0 / pi () para obtener resultados en grados.
Mike T
Sí, gracias, uso la función de grados de PostgreSQL () para eso.
mvexel
Cuchilla de carnicero. (Ni siquiera sabía que había una función de grados hasta ahora). Sería bueno concluir toda esta lógica en una llamada de función, pero tengo dificultades para conceptualizar cómo funcionaría, es decir ST_IntersectionAngle(....
Mike T
De hecho, me sorprendió que no sea una función PostGIS. Gracias por sus comentarios sobre esto.
mvexel
2

Recientemente tuve que calcular lo mismo, pero decidí un enfoque más simple y probablemente más rápido.

Para encontrar los puntos adicionales para el cálculo de acimut, solo verifico una permiría de la longitud detrás de la intersección (o después, en el raro caso de que ocurra al comienzo de la línea) usando ST_Line_Locate_Point y 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 permyriad fue arbitraria y para obtener resultados más consistentes sería mejor utilizar un desplazamiento absoluto. Para, por ejemplo, verificar 20m de antemano, cambiaría 0.0001 a 20/ST_Length(line1)y 20/ST_Length(line2)respectivamente.

lynxlynxlynx
fuente