¿Cómo calcular el ángulo de intersección de dos líneas en PostGIS?

Question

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

El punto de partida para el cálculo de ángulos en PostGIS parece ser ST_Azimut - pero que toma puntos como entrada. Mi primera idea fue tomar los puntos finales de las líneas de intersección y realizar un cálculo de acimut en ellos. Eso no es suficiente, porque la mayoría de las características de la línea no son rectas, y estoy interesado en el ángulo en la intersección. Así que lo que se me ocurrió es una operación anidada que pasa por los siguientes pasos:

1. Identifique todas las intersecciones entre las dos tablas de características lineales.
2. Crear un pequeño búfer alrededor del punto de intersección
3. Identificar los puntos en los que las características de la línea se cruzan con el exterior del búfer (tomando el primer punto si hay más de uno - en realidad sólo estoy interesado en si el ángulo es cercano a 0, 90 o 180 grados).
4. Calcule ST_Azimuth para esos dos puntos.

El SQL completo es demasiado largo para publicarlo aquí, pero lo he reunido en un gist aquí si le interesa. (Por cierto, ¿hay alguna forma mejor que traspasar todos los campos bajando por las sentencias WITH)?

Los resultados no parecen correctos, así que está claro que estoy haciendo algo mal:

EDIT He rehecho los cálculos en EPSG:3785 y los resultados son un poco diferentes pero siguen sin ser correctos:

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

Answer 1

Tuve la epifanía. Es bastante mundano. Estaba omitiendo un dato esencial para que PostGIS calculara el ángulo recto.

Lo que estaba calculando era el ángulo entre sólo los dos puntos que intersecan el pequeño exterior del búfer. Para calcular el ángulo de la intersección, tengo que calcular ambos ángulos entre los dos puntos del exterior del búfer y el punto de intersección de las dos características de línea y restarlos.

He actualizado el SQL pero aquí está lo más destacado:

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

Answer 2

Hace poco tuve que calcular lo mismo, pero me decidí por un método más sencillo y probablemente más rápido.

Para encontrar los puntos extra para el cálculo del acimut, simplemente compruebo una permirada de la longitud detrás de la intersección (o después en el raro caso de que ocurra en el mismo inicio 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 permirad era arbitraria y para obtener resultados más coherentes sería mejor utilizar un desplazamiento absoluto. Para, por ejemplo, comprobar 20 m de antemano, cambiaría 0,0001 por 20/ST_Length(line1) y 20/ST_Length(line2) respectivamente.

Answer 3

La respuesta de @mvexel es muy útil. Sin embargo, ese método no conserva la dirección de la línea. Mi método es

1. Identifique todas las intersecciones entre las dos tablas de características lineales.
2. Crea un pequeño búfer alrededor del punto de intersección.
3. Utilizar el búfer para intersectar con la otra vía. Obtener el punto final de la línea intersectada como punto 1.
4. Usando el búfer para cruzarse con la autopista. Obtener el punto final de la línea intersectada como punto 2.
5. Calcular ST_Azimuth para esos dos puntos

enlace de código

Primero: obtener puntos de intersección.

with intersections as (
select ST_Intersection(a.geometry, b.geometry) as intersection,
       a.geometry as otherway_geom,
       b.geometry as motorway_geom,
       a.id as otherway_osmid,
       b.id as motorway_osmid

from otherways a, motorways b
where ST_Intersects(a.geometry,b.geometry)),

En segundo lugar, crear un amortiguador alrededor del punto de intersección.

buffers AS (SELECT intersections.intersection, 
    ST_ExteriorRing (ST_Buffer(intersections.intersection, 1))AS extring,
    ST_Buffer(intersections.intersection, 1) as buffer,
    intersections.otherway_geom,
    intersections.motorway_geom,
    intersections.otherway_osmid,
    intersections.motorway_osmid
FROM 
    intersections),

En tercer lugar, obtener el punto final de cada línea, que preservar la dirección de linestring.

points AS(

SELECT 
    ST_Intersection(buffers.buffer, buffers.otherway_geom) as otherway_intersected_lane, 
    ST_Intersection(buffers.buffer, buffers.motorway_geom) as motorway_intersected_lane,
    st_endpoint(ST_Intersection(buffers.buffer, buffers.otherway_geom)) as point1,
    st_endpoint(ST_Intersection(buffers.buffer, buffers.motorway_geom)) as point2,
    buffers.intersection,
    buffers.extring,
    buffers.buffer,
    buffers.otherway_geom,
    buffers.motorway_geom,
    buffers.otherway_osmid,
    buffers.motorway_osmid
FROM 
    buffers)

Por último, ¡calcula el ángulo!

SELECT 

st_astext(points.point1) as point1,
st_astext(points.point2) as point2,
st_astext(points.extring) as extring,
st_astext(points.otherway_geom) as otherway_geom,
st_astext(points.motorway_geom) as motorway_geom,
abs
(
    round
    (
        degrees
        (
            ST_Azimuth
            (
                points.point2,
                points.intersection
            )

            -

            ST_Azimuth
            (
                points.point1,
                points.intersection
            )           
        )::decimal % 180.0
        ,2
    )
)AS angle_Azimuth,   
points.otherway_osmid AS otherway_id, 
points.motorway_osmid AS motorway_id 
FROM points