Dividir linhas em subconjuntos sem sobreposição com base em pontos

10

Dada uma tabela com geometria de linha e um ou mais pontos que são encaixados nessa linha em uma tabela separada, eu gostaria de dividir cada linha com um ou mais pontos de interseção em cada um dos locais onde a linha cruza um ponto.

Por exemplo, há uma linha L com três pontos de interseção, A, B e C em ordem ao longo da geometria da linha. Gostaria de retornar L como quatro geometrias distintas: do início de L a A, de A a B ao longo de L, de B a C ao longo de L e de C ao final de L.

No passado, eu usei shapely para esta tarefa, que é um problema de referência linear ( http://sgillies.net/blog/1040/shapely-recipes/ ). No entanto, isso não seria viável neste caso, que possui muitos milhões de linhas e pontos. Em vez disso, estou procurando uma solução usando o PostgreSQL / PostGIS.

Observe que os pontos são restritos a estar em uma linha. Além disso, um ponto pode validamente estar no início ou no final de uma linha; nesse caso, a linha não precisa ser dividida (a menos que haja outros pontos que não sejam coincidentes com os pontos de início ou de término da mesma linha). As linhas de subconjunto precisam manter sua direção e seus atributos, mas os atributos dos recursos de ponto não importam.

alphabetasoup
fonte

Respostas:

7

A função ST_Split PostGIS é provavelmente o que você deseja.

O PostGIS 2.2+ agora suporta Multi * geometrias no ST_Split.

Para versões mais antigas do PostGIS, leia:


Para obter uma única linha dividida por vários pontos, você pode usar algo como esta função plpgsql do empacotador multiponto . Simplifiquei apenas o caso "dividir (multi) linhas com (multi) pontos" abaixo:

DROP FUNCTION IF EXISTS split_line_multipoint(input_geom geometry, blade geometry);
CREATE FUNCTION split_line_multipoint(input_geom geometry, blade geometry)
  RETURNS geometry AS
$BODY$
    -- this function is a wrapper around the function ST_Split 
    -- to allow splitting multilines with multipoints
    --
    DECLARE
        result geometry;
        simple_blade geometry;
        blade_geometry_type text := GeometryType(blade);
        geom_geometry_type text := GeometryType(input_geom);
    BEGIN
        IF blade_geometry_type NOT ILIKE 'MULTI%' THEN
            RETURN ST_Split(input_geom, blade);
        ELSIF blade_geometry_type NOT ILIKE '%POINT' THEN
            RAISE NOTICE 'Need a Point/MultiPoint blade';
            RETURN NULL;
        END IF;

        IF geom_geometry_type NOT ILIKE '%LINESTRING' THEN
            RAISE NOTICE 'Need a LineString/MultiLineString input_geom';
            RETURN NULL;
        END IF;

        result := input_geom;           
        -- Loop on all the points in the blade
        FOR simple_blade IN SELECT (ST_Dump(ST_CollectionExtract(blade, 1))).geom
        LOOP
            -- keep splitting the previous result
            result := ST_CollectionExtract(ST_Split(result, simple_blade), 2);
        END LOOP;
        RETURN result;
    END;
$BODY$
LANGUAGE plpgsql IMMUTABLE;

-- testing
SELECT ST_AsText(split_line_multipoint(geom, blade))
    FROM (
        SELECT ST_GeomFromText('Multilinestring((-3 0, 3 0),(-1 0, 1 0))') AS geom,
        ST_GeomFromText('MULTIPOINT((-0.5 0),(0.5 0))') AS blade
        --ST_GeomFromText('POINT(-0.5 0)') AS blade
    ) AS T;

Em seguida, para criar uma geometria multiponto para cortar, use ST_Collect e crie-a manualmente a partir das entradas:

SELECT ST_AsText(ST_Collect(
  ST_GeomFromText('POINT(1 2)'),
  ST_GeomFromText('POINT(-2 3)')
));

st_astext
----------
MULTIPOINT(1 2,-2 3)

Ou agregue-o a partir de uma subconsulta:

SELECT stusps,
  ST_Multi(ST_Collect(f.the_geom)) as singlegeom
FROM (SELECT stusps, (ST_Dump(the_geom)).geom As the_geom
      FROM somestatetable ) As f
GROUP BY stusps
rcoup
fonte
Tentei ST_Split para começar e fiquei surpreso quando descobri que ele não aceitava geometria multiponto. Sua função parece preencher essa lacuna, mas infelizmente está retornando NULL para o caso multiponto de exemplo. (Funciona bem no ponto (único).) No entanto, alterei IF blade_geometry_type NOT ILIKE '% LINESTRING' THEN para IF blade_geometry_type ILIKE '% LINESTRING' THEN em sua função e obtive o resultado esperado e correto de `GEOMETRYCOLLECTION '. Ainda sou razoavelmente novo no PostGIS, portanto, essa modificação é sensata?
alphabetasoup
Desculpe, deveria ter sido IF geom_geometry_type NOT ILIKE '%LINESTRING' THEN- eu editei.
Rcoup 1/09/14
1
Ah entendo. Obrigado, esta é uma ótima solução. Você deve sugerir isso como uma contribuição ao ST_Split, para que ele possa lidar com múltiplas linhas e multipontos, se isso ainda não estiver no pipeline do PostGIS.
alphabetasoup
3
ST_Splitsuporta multi * lâminas em postgis 2.2e acima postgis.net/docs/ST_Split.html
raphael
3

Atualize para o PostGIS 2.2 , onde ST_Split foi expandido para suportar a divisão por uma multilinha, um limite multiponto ou (multi) polígono.

postgis=# SELECT postgis_version(),
                  ST_AsText(ST_Split('LINESTRING(0 0, 2 0)', 'MULTIPOINT(0 0, 1 0)'));
-[ RECORD 1 ]---+------------------------------------------------------------
postgis_version | 2.2 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
st_astext       | GEOMETRYCOLLECTION(LINESTRING(1 0,2 0),LINESTRING(0 0,1 0))
Mike T
fonte
Isto é brilhante.
alphabetasoup
Isso não funciona para o meu complicado geom: gist.github.com/ideamotor/7bd7cdee15f410ce12f3aa14ebf70177
ideamotor
ele funciona com ST_Snap, ala trac.osgeo.org/postgis/ticket/2192
ideamotor
2

Não tenho a resposta completa para você, mas ST_Line_Locate_Point usa uma linha e um ponto como argumentos e retorna um número entre 0 e 1 representando a distância ao longo da linha até a posição mais próxima do ponto.

ST_Line_Substring usa uma linha e dois números, cada um entre 0 e 1, como argumentos. Os números representam as posições na linha como distâncias fracionárias. A função retorna o segmento de linha que é executado entre essas duas posições.

Ao trabalhar com essas duas funções, você poderá conseguir o que deseja fazer.

Edward
fonte
Obrigado por isso. Na verdade, resolvi esse problema usando sua técnica e a do @rcoup. Eu dei a ele a resposta aceita devido à função que deve facilitar as coisas para os outros. Se outros quiserem seguir esse caminho, criei uma tabela temporária das linhas que possuem pontos, com uma linha para cada linha e uma parada nela. Adicionei colunas para a saída de ST_Line_Locate_Point (line.geom, pt.geom) AS L e uma função de janela: rank () OVER PARTITION BY line.id ORDER BY LR). Em seguida, deixou OUTER JOIN a tabela temporária, uma, a si mesma, b, onde a.id = b.id e a.LR = b.LR + 1 (continuação)
alphabetasoup
(continuação) A junção externa permite um CASE WHEN os campos de junção são nulos; nesse caso, ST_Line_Substring do ponto ao final da linha, ou ST_Line_Substring da referência linear do primeiro ponto à referência linear do segundo ponto (com classificação mais alta). A obtenção do segmento LA [start] é executada com um segundo SELECT, escolhendo aqueles com classificação 1 e calculando a ST_Line_Substring do ST_StartPoint da linha até a referência linear do ponto de interseção. Coloque-os na mesa, lembrando-se de manter o line.id e voilà. Felicidades.
alphabetasoup
Você pode postar esta resposta como uma resposta no código, por favor? Eu gostaria de olhar para essa opção e também sou um pouco novato no SQL.
Phil Donovan
1
@PhilDonovan: done.
alphabetasoup
2

Me pediram isso duas vezes agora, desculpe pelo atraso. É improvável que seja considerada uma solução concisa; Eu o escrevi um pouco mais abaixo na curva de aprendizado do que atualmente. Todas as dicas são bem-vindas, até as estilísticas.

--Inputs:
--walkingNetwork = Line features representing edges pedestrians can walk on
--stops = Bus stops
--NOTE: stops.geom is already constrained to be coincident with line features
--from walkingNetwork. They may be on a vertex or between two vertices.

--This series of queries returns a version of walkingNetwork, with edges split
--into separate features where they intersect stops.

CREATE TABLE tmp_lineswithstops AS (
    WITH subq AS (
        SELECT
        ST_Line_Locate_Point(
            roads.geom,
            ST_ClosestPoint(roads.geom, stops.geom)
        ) AS LR,
        rank() OVER (
            PARTITION BY roads.gid
            ORDER BY ST_Line_Locate_Point(
                roads.geom,
                ST_ClosestPoint(roads.geom, stops.geom)
            )
        ) AS LRRank,
        ST_ClosestPoint(roads.geom, stops.geom),
        roads.*
        FROM walkingNetwork AS roads
        LEFT OUTER JOIN stops
        ON ST_Distance(roads.geom, stops.geom) < 0.01
        WHERE ST_Equals(ST_StartPoint(roads.geom), stops.geom) IS false
        AND ST_Equals(ST_EndPoint(roads.geom), stops.geom) IS false
        ORDER BY gid, LRRank
    )
    SELECT * FROM subq
);

-- Calculate the interior edges with a join
--If the match is null, calculate the line to the end
CREATE TABLE tmp_testsplit AS (
    SELECT
    l1.gid,
    l1.geom,
    l1.lr AS LR1,
    l1.st_closestpoint AS LR1geom,
    l1.lrrank AS lr1rank,
    l2.lr AS LR2,
    l2.st_closestpoint AS LR2geom,
    l2.lrrank AS lr2rank,
    CASE WHEN l2.lrrank IS NULL -- When the point is the last along the line
        THEN ST_Line_Substring(l1.geom, l1.lr, 1) --get the substring line to the end
        ELSE ST_Line_Substring(l1.geom, l1.lr, l2.lr) --get the substring between the two points
    END AS sublinegeom
    FROM tmp_lineswithstops AS l1
    LEFT OUTER JOIN tmp_lineswithstops AS l2
    ON l1.gid = l2.gid
    AND l2.lrrank = (l1.lrrank + 1)
);

--Calculate the start to first stop edge
INSERT INTO tmp_testsplit (gid, geom, lr1, lr1geom, lr1rank, lr2, lr2geom, lr2rank, sublinegeom)
SELECT gid, geom,
0 as lr1,
ST_StartPoint(geom) as lr1geom,
0 as lr1rank,
lr AS lr2,
st_closestpoint AS lr2geom,
lrrank AS lr2rank,
ST_Line_Substring(l1.geom, 0, lr) AS sublinegeom --Start to point
FROM tmp_lineswithstops AS l1
WHERE l1.lrrank = 1;

--Now match back to the original road features, both modified and unmodified
CREATE TABLE walkingNetwork_split AS (
    SELECT
    roadssplit.sublinegeom,
    roadssplit.gid AS sgid, --split-gid
    roads.*
    FROM tmp_testsplit AS roadssplit
    JOIN walkingNetwork AS r
    ON r.gid = roadssplit.gid
    RIGHT OUTER JOIN walkingNetwork AS roads --Original edges with null if unchanged, original edges with split geom otherwise
    ON roads.gid = roadssplit.gid
);

--Now update the necessary columns, and drop the temporary columns
--You'll probably need to work on your own length and cost functions
--Here I assume it's valid to just multiply the old cost by the fraction of
--the length the now-split line represents of the non-split line
UPDATE walkingNetwork_split
SET geom = sublinegeom,
lengthz = lengthz*(ST_Length(sublinegeom)/ST_Length(geom)),
walk_seconds_ft = walk_seconds_ft*(ST_Length(sublinegeom)/ST_Length(geom)),
walk_seconds_tf = walk_seconds_tf*(ST_Length(sublinegeom)/ST_Length(geom))
WHERE sublinegeom IS NOT NULL
AND ST_Length(sublinegeom) > 0;
ALTER TABLE walkingNetwork_split
DROP COLUMN sublinegeom,
DROP COLUMN sgid;

--Drop intermediate tables
--You probably could use actual temporary tables;
--I prefer to have a sanity check at each stage
DROP TABLE IF EXISTS tmp_testsplit;
DROP TABLE IF EXISTS tmp_lineswithstops;

--Assign the edges a new unique id, so we can use this as source/target columns in pgRouting
ALTER TABLE walkingNetwork_split
DROP COLUMN IF EXISTS fid;
ALTER TABLE walkingNetwork_split
ADD COLUMN fid INTEGER;
CREATE SEQUENCE roads_seq;
UPDATE walkingNetwork_split
SET fid = nextval('roads_seq');
ALTER TABLE walkingNetwork_split
ADD PRIMARY KEY ("fid");
alphabetasoup
fonte
0

Eu quero expandir as respostas acima da perspectiva de um iniciante. Nesse cenário, você tem uma série de pontos e observa como usá-los como uma "lâmina" para cortar linhas em segmentos. Este exemplo inteiro pressupõe que você primeiro inseriu seus pontos na linha e que os pontos têm o atributo de ID exclusivo da sua linha ajustada. Eu uso 'column_id "para representar o ID exclusivo da linha.

Primeiro , você deseja agrupar seus pontos em multipontos quando mais de uma lâmina cair em uma linha. Caso contrário, a função split_line_multipoint atua como a função ST_Split, que não é o resultado desejado.

CREATE TABLE multple_terminal_lines AS
SELECT ST_Multi(ST_Union(the_geom)) as the_geom, a.matched_alid
FROM    point_table a
        INNER JOIN
        (
            SELECT  column_id
            FROM    point_table
            GROUP   BY column_id
            HAVING  COUNT(*) > 1
        ) b ON a.column_id = b.column_id
GROUP BY a.column_id;

Então , você deseja dividir sua rede com base nesses multipontos.

CREATE TABLE split_multi AS
SELECT (ST_Dump(split_line_multipoint(ST_Snap(a.the_geometry, b.the_geom, 0.00001),b.the_geom))).geom as the_geom
FROM line_table a
JOIN multple_terminal_lines b 
ON a.column_id = b.column_id;


Repita as etapas 1 e 2 com suas linhas que possuem apenas um ponto de interseção. Para fazer isso, você deve atualizar o código da etapa 1 para 'TENDO CONTAGEM (*) = 1'. Renomeie as tabelas de acordo.


Em seguida , crie uma tabela de linhas duplicada e exclua as entradas com pontos nelas.

CREATE TABLE line_dup AS
SELECT * FROM line_table;
-- Delete shared entries
DELETE FROM line_dup
WHERE column_id in (SELECT DISTINCT column_id FROM split_single) OR column_id in (SELECT DISTINCT column_id FROM split_multi) ;


Por fim , junte suas três tabelas usando UNION ALL:

CREATE TABLE line_filtered AS 
SELECT the_geom
FROM split_single
UNION ALL 
SELECT the_geom
FROM split_multi
UNION ALL 
SELECT the_geom
FROM line_dup;

BAM!

Mtap1
fonte