¿Identificar puntos en el lado izquierdo o derecho de la calle utilizando PostGIS?

Question

El escenario:

Tengo puntos (direcciones) y cadenas multilíneas (calles), es decir ap y street tablas en mi db PostgreSQL 9.5. Tabla ap contiene 3 filas (puntos de muestra) mientras que street contiene 85 filas. El escenario de ejemplo se muestra a continuación.

Lo que necesito:

Para cada punto ap Necesito identificar si está en el "lado izquierdo" o en el "lado derecho" del street .

El código:

En GIS SE, ya he encontrado algunas preguntas relevantes. Por ejemplo:

• ¿Detectar si un punto está en el lado izquierdo o derecho de una línea en PostGIS? y
• ¿Determinar en qué lado se encuentra un punto de la línea utilizando PostGIS?

El siguiente código es una adaptación del primero enlace .

With cp as
(
Select
    a.gid ap_id,
    a.geom ap,
    ST_SetSRID(ST_LineMerge(b.geom), 3044) line,
    ST_SetSRID(ST_ClosestPoint(b.geom, a.geom), 3044) p
from
    ap a
Left Join street b on
ST_Dwithin(ST_SetSRID(b.geom, 3044), ST_SetSRID(a.geom, 3044), 25)
Order by a.gid, ST_Distance(b.geom, a.geom)
)
, h as
(
Select
    ap_id, ap,
    line, p,
        ST_MakeLine(p, ap) vec,
        ST_MakeLine(p,
        ST_LineInterpolatePoint(line,
            ST_LineLocatePoint(line, p) * 1.01)) seg
From cp
Order by ap_id
)
Select
    ap_id,
    degrees(ST_Azimuth(ST_StartPoint(vec), ST_EndPoint(vec))) 
        - degrees(ST_Azimuth(ST_StartPoint(seg), ST_EndPoint(seg))) diff
from h
order by ap_id

La salida:

El código anterior (para el escenario de ejemplo) funciona y la salida es:

ap_id    diff
1       -269.999999988824
2        89.9999999946886
3        269.999999980956

Error:

Sin embargo, el código para otro conjunto de datos de puntos (87 filas), arroja el siguiente error:

ERROR:  line_interpolate_point: 2nd arg isn't within [0,1]

Además, el código anterior no parece funcionar para los casos en los que la calle es de este a oeste y tiene un rumbo de 90 grados (Referencia, aquí ). Necesito realizar esta tarea de identificación de puntos para un gran conjunto de datos. Por lo tanto, estoy buscando sugerencias para superar este problema. ¿Hay alguna forma alternativa de identificar si los puntos están en el lado izquierdo o en el lado derecho de la calle?

Answer 1

Esto hace uso de su consulta actual, sin embargo he puesto ST_ShorttestLine del ST_ClosestPoint y han añadido ST_LineCrossingDirection a la salida. Yo esperaría que devolviera 1 o -1 para la izquierda o la derecha. Sin embargo, es posible que necesite que se extienda la línea más corta para asegurarse de que se cruza. En geometrías muy simples esto parece funcionar a menos que el punto esté por delante de la línea (por ejemplo, el punto más cercano es el último punto de la calle)

With cp as
(
Select
    a.gid ap_id,
    a.geom ap,
    ST_SetSRID(ST_LineMerge(b.geom), 3044) line,
    ST_SetSRID(ST_ShortestLine(b.geom, a.geom), 3044) vec
from
    ap a
Left Join street b on
ST_Dwithin(ST_SetSRID(b.geom, 3044), ST_SetSRID(a.geom, 3044), 25)
Order by a.gid, ST_Distance(b.geom, a.geom)
)
Select *,
    ST_LineCrossingDirection(vec,line) side
from cp;

Teniendo en cuenta que es probable que el error de punto flotante comience a causar algunos problemas y para tratar de abordar los puntos principales, he ajustado la consulta para extender las líneas añadiendo puntos. Las extensiones se basan en esto responder . Esto parece afectar a los puntos que están en el eje, pero puede que no sean importantes para usted.

He dejado mis datos de prueba en la consulta.

WITH ap AS (
SELECT gid
    ,ST_SetSRID(ST_MakePoint(x,y),3044) geom
    ,description
FROM (VALUES
     (1, -1, 0, 'Left Before')
    ,(2, 0, 1, 'Left Middle')
    ,(3, 2, 3, 'Left After')
    ,(4, 0, -1, 'Right Before')
    ,(5, 1, 0, 'Right Middle')
    ,(6, 3, 2, 'Right After')
    ,(7, -1, -1, 'On Before')
    ,(8, 1, 1, 'On Middle')
    ,(9, 3, 3, 'On After')
    ) ap(gid, x, y, description)
),
street AS (
SELECT ST_SetSRID(ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(2,2)),3044) geom
),
--With 
       cp as
(
Select
    a.gid ap_id,
    a.geom ap,
    a.description,
    ST_SetSRID( -- Line with an extension added to it
        ST_AddPoint(
            ST_MakeLine(
                a.geom,
                ST_ClosestPoint(b.geom,a.geom)
            ),
            ST_Translate(
                a.geom, 
                sin(ST_Azimuth(a.geom, ST_ClosestPoint(b.geom,a.geom))) * ST_Distance(a.geom, ST_ClosestPoint(b.geom,a.geom)) * 1.01,
                cos(ST_Azimuth(a.geom, ST_ClosestPoint(b.geom,a.geom))) * ST_Distance(a.geom, ST_ClosestPoint(b.geom,a.geom)) * 1.01
            )
        )
    , 3044) vec,
    ST_SetSRID(
        ST_AddPoint(
            ST_LineMerge(b.geom),
            ST_Translate(
                ST_PointN(ST_LineMerge(b.geom),-2),
                sin(ST_Azimuth(ST_PointN(ST_LineMerge(b.geom),-2), ST_PointN(ST_LineMerge(b.geom),-1))) * (ST_Distance(ST_PointN(ST_LineMerge(b.geom),-1), ST_PointN(ST_LineMerge(b.geom),-2)) *1.1),
                cos(ST_Azimuth(ST_PointN(ST_LineMerge(b.geom),-2), ST_PointN(ST_LineMerge(b.geom),-1))) * (ST_Distance(ST_PointN(ST_LineMerge(b.geom),-1), ST_PointN(ST_LineMerge(b.geom),-2)) *1.1)
            )
        )
    , 3044) line
from
    ap a
Left Join street b on
ST_Dwithin(ST_SetSRID(b.geom, 3044), ST_SetSRID(a.geom, 3044), 25)
Order by a.gid, ST_Distance(b.geom, a.geom)
)
Select *,
    ST_LineCrossingDirection(line,vec) side
from cp;

Dando el siguiente resultado

Answer 2

Comprobar si el punto está en el lado izquierdo o derecho con respecto a la dirección de la calle. Basta con crear un buffer de línea izquierda o derecha con la misma distancia al punto. Si el punto se cruza con el buffer entonces debe estar a la izquierda, sino a la derecha. Este es el código de ejemplo. Puedes reemplazar intersecciones con cualquier otra función.

SELECT  
 CASE
   WHEN ST_Intersects (
      ST_Buffer(line_geom, 
         ST_Distance(line_geom, point_geom)*1.1, 'side=left endcap=flat'),
      point_geom) = TRUE
  THEN 'T'
  ELSE 'F'
  END AS isleft
FROM some_table;