Fractionner les lignes en sous-ensembles non superposés en fonction des points

10

Étant donné une table avec une géométrie de ligne et un ou plusieurs points qui sont alignés sur cette ligne dans une table distincte, je voudrais diviser chaque ligne avec un ou plusieurs points d'intersection à chacun des emplacements où la ligne coupe un point.

Par exemple, il existe une ligne, L, avec trois points d'intersection, A, B et C dans l'ordre le long de la géométrie de la ligne. Je voudrais renvoyer L comme quatre géométries distinctes: du début de L à A, de A à B le long de L, de B à C le long de L et de C à la fin de L.

Dans le passé, j'ai utilisé galbée pour cette tâche, qui est un problème de référencement linéaire ( http://sgillies.net/blog/1040/shapely-recipes/ ). Cependant, cela ne serait pas réalisable dans ce cas, qui a plusieurs millions de lignes et de points. Au lieu de cela, je recherche une solution utilisant PostgreSQL / PostGIS.

Notez que les points sont contraints d'être sur une ligne. De plus, un point peut être valablement au début ou à la fin d'une ligne, auquel cas la ligne n'a pas besoin d'être divisée (sauf s'il y a d'autres points qui ne coïncident pas avec les points de début ou de fin de la même ligne). Les lignes de sous-ensemble doivent conserver leur direction et leurs attributs, mais les attributs des entités ponctuelles n'ont pas d'importance.

alphabetasoup
la source

Réponses:

7

La fonction ST_Split PostGIS est probablement ce que vous voulez.

PostGIS 2.2+ prend désormais en charge les géométries Multi * dans ST_Split.

Pour les anciennes versions de PostGIS, lisez la suite:


Pour obtenir une seule ligne divisée en plusieurs points, vous pouvez utiliser quelque chose comme cette fonction plpgsql d' encapsuleur multipoint . Je l'ai simplifié juste au cas "lignes (multi) divisées avec (multi) points" ci-dessous:

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;

Ensuite, pour créer une géométrie multipoint à découper, utilisez ST_Collect et créez-la manuellement à partir des entrées:

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

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

Ou agrégez-le à partir d'une sous-requête:

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
la source
J'ai essayé ST_Split pour commencer, et j'ai été surpris quand j'ai découvert qu'il n'acceptait pas la géométrie multipoint. Votre fonction semble combler cette lacune, mais malheureusement elle retourne NULL pour l'exemple de cas multipoint. (Cela fonctionne bien sur le point (unique).) Cependant, j'ai changé IF blade_geometry_type NOT ILIKE '% LINESTRING' THEN à IF blade_geometry_type ILIKE '% LINESTRING' THEN dans votre fonction et j'ai obtenu le résultat attendu et correct de `GEOMETRYCOLLECTION '. Cependant, je suis encore relativement nouveau sur PostGIS, cette modification est-elle donc raisonnable?
alphabetasoup
Désolé, aurait dû être IF geom_geometry_type NOT ILIKE '%LINESTRING' THEN- je l'ai édité.
rcoup
1
Ah, je vois. Merci, c'est une excellente solution. Vous devez suggérer cela comme une contribution à ST_Split afin qu'il puisse gérer les multilignes et les multipoints, si ce n'est pas déjà dans le pipeline PostGIS.
alphabetasoup
3
ST_Splitprend en charge les lames multi * dans postgis 2.2et au-dessus de postgis.net/docs/ST_Split.html
raphael
3

Mise à niveau vers PostGIS 2.2 , où ST_Split a été développé pour prendre en charge la division par une frontière multiligne, multipoint ou (multi) polygone.

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
la source
C'est génial.
alphabetasoup
Cela ne fonctionne pas pour mon geom compliqué: gist.github.com/ideamotor/7bd7cdee15f410ce12f3aa14ebf70177
ideamotor
cela fonctionne avec ST_Snap, ala trac.osgeo.org/postgis/ticket/2192
ideamotor
2

Je n'ai pas toute la réponse pour vous, mais ST_Line_Locate_Point prend une ligne et un point comme arguments, et renvoie un nombre compris entre 0 et 1 représentant la distance le long de la ligne jusqu'à la position la plus proche du point.

ST_Line_Substring prend une ligne et deux nombres, chacun entre 0 et 1, comme arguments. Les nombres représentent les positions sur la ligne sous forme de distances fractionnaires. La fonction renvoie le segment de ligne qui s'étend entre ces deux positions.

En travaillant avec ces deux fonctions, vous devriez pouvoir réaliser ce que vous voulez faire.

Edward
la source
Merci pour cela. J'ai en fait résolu ce problème en utilisant votre technique ainsi que celle de @rcoup. Je lui ai donné la réponse acceptée en raison de la fonction qui devrait faciliter la tâche des autres. Si d'autres veulent suivre ce chemin, j'ai créé un tableau temporaire des lignes qui ont des points sur eux, avec une ligne pour chaque ligne et un arrêt qui s'y trouve. J'ai ajouté des colonnes pour la sortie de ST_Line_Locate_Point (line.geom, pt.geom) AS L et une fonction de fenêtre: rank () OVER PARTITION BY line.id ORDER BY LR). Puis LEFT OUTER JOIN la table temporaire, a, à elle-même, b, où a.id = b.id et a.LR = b.LR + 1 (suite)
alphabetasoup
(suite) La jointure externe permet un CASE QUAND les champs de jointure sont nuls, auquel cas ST_Line_Substring du point à la fin de la ligne, sinon ST_Line_Substring de la référence linéaire du premier point à la référence linéaire du deuxième point (avec un rang plus élevé). L'obtention du segment [start] LA est ensuite effectuée avec un second SELECT, il suffit de sélectionner ceux avec un rang de 1 et de calculer la ST_Line_Substring du ST_StartPoint de la ligne jusqu'à la référence linéaire du point d'intersection. Mettez-les dans le tableau, en vous souvenant de garder le line.id, et voilà. À votre santé.
alphabetasoup
Pouvez-vous poster cette réponse en tant que réponse en code s'il vous plaît? Je voudrais examiner cette option ainsi que je suis un peu novice en SQL.
Phil Donovan
1
@PhilDonovan: terminé.
alphabetasoup
2

On me l'a demandé deux fois maintenant, donc désolé pour le retard. Il est peu probable que cela soit considéré comme une solution laconique; Je l'ai écrit un peu plus loin dans la courbe d'apprentissage que je ne le suis actuellement. Tous les conseils sont les bienvenus, même ceux stylistiques.

--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
la source
0

Je veux développer les réponses ci-dessus du point de vue d'un débutant. Dans ce scénario, vous avez une série de points et vous regardez pour les utiliser comme une "lame" pour couper des lignes en segments. Cet exemple suppose que vous avez d'abord accroché vos points à la ligne et que les points ont l'attribut ID unique de leur ligne accrochée. J'utilise 'column_id "pour représenter l'ID unique de la ligne.

Tout d'abord , vous souhaitez regrouper vos points en multipoints lorsque plus d'une lame tombe sur une ligne. Sinon, la fonction split_line_multipoint agit comme la fonction ST_Split, qui n'est pas le résultat souhaité.

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;

Ensuite , vous souhaitez diviser votre réseau en fonction de ces multipoints.

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;


Répétez les étapes 1 et 2 avec vos lignes qui n'ont qu'un seul point d'intersection. Pour ce faire, vous devez mettre à jour le code de l'étape 1 en 'HAVING COUNT (*) = 1'. Renommez les tables en conséquence.


Ensuite , créez un tableau de lignes en double et supprimez les entrées avec des points dessus.

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


Enfin , réunissez vos trois tables en utilisant 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
la source