In questo blog post vedremo come effettuare la differenza tra linee sovrapposte, lo faremo con due strumenti: QGIS e SpatiaLite.

Dalla guida online di PostGIS:

Ripasso:

Dalla guida sembra facile effettuare la differenza; vediamo un esempio pratico: in QGIS creo due linee (due semplici segmenti sovrapposti parzialmente) in due layer differenti:

La differenza tra linea A e linea B che mi aspetto è:

da Processing di QGIS cerco l’algoritmo Differenza e imposto la maschera di input:

il risultato che ottengo è palesemente ERRATO:

ottengo nuovamente la linea A, PERCHÈ?

Leggendo la definizione di differenza nella guida di PostGIS troviamo:

Si può pensare a questo come GeometryA - ST_Intersection (A, B)

Faccio l’intersezione tra linea A e linea B in QGIS, ottengo anche questa volta un risultato (apparentemente) errato:

cioè, non restituisce nulla, insieme VUOTO! (PS. le due linee sono state create usando gli snap, quindi la sovrapposizione c’è, ma NON è sufficiente per le operazioni di overlapping).

In sintesi: per poter ottenere risultati corretti occorre che le linee abbiano in corrispondenza della sovrapposizione:

• dei nodi di disegno
• e i nodi devono essere perfettamente snappati agli estremi delle linee sovrapposte.

Come procedere

in QGIS:

Usare (con attenzione) algoritmo Aggancia geometria al vettore usando l’opzione corretta. Questo algoritmo aggiungerà dei nodi, nel nostro caso quelli che ci servono.

oppure tramite il field calc, in aggiornamento della geometria, usare l’espressione:

 make_line(  
 start_point($geometry), 
 end_point(geometry(get_feature_by_id('linea_a',1))),
 end_point($geometry))
 )

dopo aver aggiunto i due nodi di disegno, alla linea a e linea b, occorre snappare e quindi riutilizzare nuovamente l’algoritmo Aggancia geometria al vettore;

infine usare l’algoritmo Differenza e otterrete, finalmente, il tratto di geometria desiderato.

in SpatiaLite e SQL:

Proiettare il vertice terminale della linea a sulla linea b con la seguente query:

DROP TABLE IF EXISTS "points_over_lines_linea_a";

CREATE TABLE "points_over_lines_linea_a" AS
SELECT ROW_NUMBER() OVER( ORDER BY "a.pkuid" ) AS pk, 
ST_ClosestPoint(c.geom, st_translate(a.geom,0.001,0.001,0)) AS geom
FROM (SELECT StartPoint(geom) as geom FROM "linea_a") a, "linea_b" c;

SELECT RecoverGeometryColumn('points_over_lines_linea_a','geom',32632,'POINT','XY');

SELECT CreateSpatialIndex('points_over_lines_linea_a', 'geom');

ottenendo questo:

DROP TABLE IF EXISTS "points_over_lines_linea_b";

CREATE TABLE "points_over_lines_linea_b" AS
SELECT ROW_NUMBER() OVER( ORDER BY "a.pkuid" ) AS pk, 
ST_ClosestPoint(c.geom, st_translate(a.geom,0.001,0.001,0)) AS geom
FROM (SELECT StartPoint(geom) as geom FROM "linea_a") a, "linea_b" c;

SELECT RecoverGeometryColumn('points_over_lines_linea_b','thegeom',32632,'POINT','XY');

SELECT CreateSpatialIndex('points_over_lines_linea_b', 'thegeom');

ottenendo:

Questi due punti devono essere aggiungi alle geometrie delle linee (linea a e linea b)

-- aggiorna geotabella 'linea_b' aggiunge nodi alla geometria
UPDATE "linea_b" SET geom=
    CastToLinestring(CastToMulti(
    RemoveRepeatedPoints(
    ST_Snap( "linea_b".geom,
    (SELECT CastToMultipoint(st_collect(b.geom)) 
    FROM "points_over_lines_linea_b" as b
    WHERE b.pk = "linea_b".pkuid
    GROUP BY b.pk) , 0.01 ), 0.01 )
 )) 
WHERE EXISTS(
            SELECT 1 FROM "points_over_lines_linea_b" as b
            WHERE  b.pk = "linea_b".pkuid limit 1);

-- aggiorna geotabella 'linea_a' aggiunge nodi alla geometria
UPDATE "linea_a" SET geom=
    CastToLinestring(CastToMulti(
    RemoveRepeatedPoints(
    ST_Snap( "linea_a".geom,
    (SELECT CastToMultipoint(st_collect(b.geom)) 
    FROM "points_over_lines_linea_a" as b
    WHERE b.pk = "linea_a".pkuid
    GROUP BY b.pk) , 0.01 ), 0.01 )
 )) 
WHERE EXISTS(
            SELECT 1 FROM "points_over_lines_linea_a" as b
            WHERE  b.pk = "linea_a".pkuid limit 1);

Ora siamo pronti per fare la differenza tra le due linee:

CREATE TABLE "diff_a_b" AS
SELECT St_Difference(a.geom, b.geom) as geom
FROM linea_a a, linea_b b;

SELECT RecoverGeometryColumn('diff_a_b','geom',32632,'LINESTRING','XY');

NOTA FINALE: Cinque sono i gradi per giungere alla saggezza: tacere, ascoltare, ricordare, agire, studiare.
(Proverbio arabo)

Database e progetto QGIS

Riferimenti:

• QGIS : https://qgis.org/it/site/
• SpatiaLite : https://www.gaia-gis.it/fossil/libspatialite/index
• PostGIS : https://postgis.net/
• SQL : https://it.wikipedia.org/wiki/Structured_Query_Language

Se il blog post vi è piaciuto cliccate su ‘Mi piace’, grazie!!!
if you liked the blog post click on ‘Like’, thank you !!!

Se il post vi è stato utile contribuite a mantenerlo aggiornato PayPal