Podziel linie na nienakładające się podzbiory na podstawie punktów

10

Biorąc pod uwagę tabelę z geometrią linii i jeden lub więcej punktów, które są przyciągane do tej linii w osobnej tabeli, chciałbym podzielić każdą linię jednym lub więcej przecinającymi się punktami w każdym z miejsc, w których linia przecina punkt.

Na przykład istnieje linia L z trzema przecinającymi się punktami A, B i C w kolejności wzdłuż geometrii linii. Chciałbym zwrócić L jako cztery różne geometrie: od początku L do A, od A do B wzdłuż L, od B do C wzdłuż L i od C do końca L.

W przeszłości używałem zgrabnie do tego zadania, które jest liniowym problemem odniesienia ( http://sgillies.net/blog/1040/shapely-recipes/ ). Nie byłoby to jednak możliwe w tym przypadku, który ma wiele milionów linii i punktów. Zamiast tego szukam rozwiązania wykorzystującego PostgreSQL / PostGIS.

Zauważ, że punkty są ograniczone do linii. Ponadto punkt może znajdować się na początku lub na końcu linii, w którym to przypadku linia nie musi być dzielona (chyba że istnieją inne punkty, które nie są zbieżne z punktami początkowymi lub końcowymi tej samej linii). Linie podzestawu muszą zachować swój kierunek i atrybuty, ale atrybuty cech punktowych nie mają znaczenia.

zupa alfabetyczna
źródło

Odpowiedzi:

7

Funkcja ST_Split PostGIS jest prawdopodobnie tym, czego chcesz.

PostGIS 2.2+ obsługuje teraz geometrie Multi * w ST_Split.

W przypadku starszych wersji PostGIS czytaj dalej:


Aby uzyskać podział jednej linii przez wiele punktów, możesz użyć czegoś takiego jak ta funkcja wielopunktowego opakowania plpgsql. Uprościłem to tylko do poniższego przypadku „dzielenia (wielu) linii z (wieloma) punktami”:

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;

Następnie, aby utworzyć wielopunktową geometrię do wycięcia, użyj ST_Collect i albo utwórz ją ręcznie z danych wejściowych:

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

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

Lub agreguj go z podzapytania:

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
źródło
Próbowałem na początku ST_Split i byłem zaskoczony, gdy stwierdziłem, że nie akceptuje geometrii wielopunktowej. Twoja funkcja wydaje się wypełniać tę lukę, ale niestety zwraca wartość NULL dla przykładowego przypadku wielopunktowego. (Działa dobrze na (pojedynczy) pkt.), Jednak zmieniłem IF blade_geometry_type NIE ILIKE „% LineString”, a następnie do IF blade_geometry_type iLike „% LINESTRING”, a następnie w swojej funkcji, ale i oczekiwanego prawidłowego `GeometryCollection” wynik. Nadal jestem jednak dość nowy w PostGIS, więc czy ta modyfikacja jest sensowna?
alphabetasoup
Przepraszam, powinienem był IF geom_geometry_type NOT ILIKE '%LINESTRING' THEN- ja to zredagowałem.
rcoup
1
O, rozumiem. Dzięki, to świetne rozwiązanie. Powinieneś zasugerować to jako wkład do ST_Split, aby mógł obsługiwać multilinię i wielopunkt, jeśli nie jest to jeszcze w potoku PostGIS.
alphabetasoup
3
ST_Splitobsługuje multi * ostrzy postgis 2.2i powyżej postgis.net/docs/ST_Split.html
Raphael
3

Uaktualnij do PostGIS 2.2 , gdzie ST_Split został rozszerzony, aby obsługiwać podział przez granicę multilinii, wielopunktów lub (wielu) wielokątów.

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.
źródło
To jest genialne.
alphabetasoup
To nie działa dla mojej skomplikowanej geom: gist.github.com/ideamotor/7bd7cdee15f410ce12f3aa14ebf70177
ideamotor
współpracuje z ST_Snap, ala trac.osgeo.org/postgis/ticket/2192
ideamotor
2

Nie mam dla ciebie całej odpowiedzi, ale ST_Line_Locate_Point przyjmuje argumenty i linię i zwraca liczbę od 0 do 1 reprezentującą odległość wzdłuż linii do pozycji najbliższej punktu.

ST_Line_Substring przyjmuje jako argument linię i dwie liczby, każda od 0 do 1. Liczby reprezentują pozycje na linii jako odległości ułamkowe. Funkcja zwraca segment linii biegnący między tymi dwiema pozycjami.

Pracując z tymi dwiema funkcjami, powinieneś być w stanie osiągnąć to, co chcesz zrobić.

Edward
źródło
Dzięki za to. Naprawdę rozwiązałem ten problem za pomocą twojej techniki, jak również z @rcoup. Udzieliłem mu akceptowanej odpowiedzi ze względu na funkcję, która powinna ułatwić innym. Jeśli inni chcą zejść tą ścieżką, stworzyłem tymczasową tabelę linii, które mają na nich punkty, z wierszem dla każdej linii i jednym przystankiem na niej. Dodałem kolumny dla danych wyjściowych ST_Line_Locate_Point (line.geom, pt.geom) AS L i funkcji okna: rank () OVER PARTITION BY line.id ORDER BY LR). Następnie LEWA ZEWNĘTRZNA DOŁĄCZ do tabeli tymczasowej, a, do siebie, b, gdzie a.id = b.id i a.LR = b.LR + 1 (ciąg dalszy)
alfabetasoup
(ciąg dalszy) Sprzężenie zewnętrzne pozwala na PRZYPADEK, GDY pola łączenia są puste, w takim przypadku ST_Line_Substring od punktu do końca linii, w przeciwnym razie ST_Line_Substring od liniowego odniesienia pierwszego punktu do liniowego odniesienia drugiego punktu (z wyższą rangą). Uzyskanie [początkowego] odcinka LA jest następnie wykonywane za pomocą drugiego WYBORU, po prostu wybierając te z rangą 1 i obliczając ST_Line_Substring z ST_StartPoint linii do odniesienia liniowego punktu przecięcia. Włóż je do stołu, pamiętając o zachowaniu line.id i voilà. Twoje zdrowie.
alphabetasoup
Czy możesz zamieścić tę odpowiedź jako odpowiedź w kodzie? Chciałbym spojrzeć na tę opcję, a także jestem trochę nowicjuszem w SQL.
Phil Donovan,
1
@PhilDonovan: gotowe.
alphabetasoup
2

Zostałem już o to poproszony dwa razy, więc przepraszam za opóźnienie. Jest mało prawdopodobne, aby uznać to za zwięzłe rozwiązanie; Napisałem to, gdy jestem nieco dalej niż w trakcie nauki. Wszelkie wskazówki są mile widziane, nawet te stylistyczne.

--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");
zupa alfabetyczna
źródło
0

Chcę rozwinąć powyższe odpowiedzi z perspektywy początkującego. W tym scenariuszu masz serię punktów i patrzysz, aby użyć ich jako „ostrza” do cięcia linii na segmenty. Cały ten przykład zakłada, że ​​najpierw przyciągnąłeś punkty do linii i że punkty mają unikalny atrybut ID z ich linii przyciągania. Używam „column_id” do reprezentowania unikalnego identyfikatora linii.

Po pierwsze , chcesz grupie swoje punkty na multipoints gdy więcej niż jedno ostrze spada na linię. W przeciwnym razie funkcja split_line_multipoint działa jak funkcja ST_Split, co nie jest pożądanym rezultatem.

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;

Następnie chcesz podzielić sieć na podstawie tych wielopunktów.

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;


Powtórz kroki 1 i 2 dla linii, które mają tylko jeden punkt przecięcia. Aby to zrobić, należy zaktualizować kod z kroku 1 do „HAVING COUNT (*) = 1”. Zmień nazwy tabel odpowiednio.


Następnie utwórz zduplikowaną tabelę linii i usuń wpisy z punktami na nich.

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


Na koniec połącz trzy tabele razem, używając 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
źródło