Retrieving the complementary difference of 2 geometries

[mj10777]

(Not sure if 'complementary difference' is the right terminology)

For a historical street Database the attempt is being made to adapt the street-segments (portions of the street mostly between 2 intersections)

- so that those portions crossing border-zone can be split

-- into the portions still accessible to the public (from either side, East or West)

-- within the border zone will be inactive (valid_until 1961-03-12, valid_since 1990-07-01)

A Polygon (called 'soldner_polygon') exists displaying the Border Zone (Outer, Inner Wall). [Blue area]

The yellow dotted line show such a segment (called 'soldner_segment'), which only contains 2 points (Start/End).

The goal is to split the yellow line into 3 portions.

With

CastToMultiLineString(ST_Intersection(b_source.soldner_polygon, a_source.soldner_segment)) AS soldner_segment_mauer

the portion inside the Polygon can be extracted (soldner_segment_mauer).

What I have not figured out is how to:

- extract from 'soldner_segment' the portion not contained in 'soldner_segment_mauer'

My first attempt was:

soldner_segment_diff=CastToMultiLineString(ST_Difference(CastToMultiLineString(soldner_segment),soldner_segment_mauer))

the returned result was the same line, but also containing the Start/End points of 'soldner_segment_mauer' together with the original Start/End points.

No doubt the solution has been staring at my without blinking (and thus not noticed)...

Will try ST_SymDifference as the next attempt tomorrow.

Mark

[a.fu...@lqt.it]

Hi Mark,

I wasn't able to exactly reproduce your problem;
Something doesn't add up.

I started by creating several tables attempting
to closely mimic your own names as far as possible:

CREATE TABLE mauer (id INTEGER PRIMARY KEY);
SELECT AddGeometryColumn('mauer', 'geom', -1,
'POLYGON', 'XY');
CREATE TABLE segment (id INTEGER PRIMARY KEY);
SELECT AddGeometryColumn('segment', 'geom', -1,
'LINESTRING', 'XY');
CREATE TABLE segment_mauer (id INTEGER PRIMARY KEY);
SELECT AddGeometryColumn('segment_mauer', 'geom', -1,
'LINESTRING', 'XY');
CREATE TABLE segment_no_mauer (id INTEGER PRIMARY KEY);
SELECT AddGeometryColumn('segment_no_mauer', 'geom', -1,
'MULTILINESTRING', 'XY');

Then I've inserted a rectangle (" mauer") and a simple
segment:

INSERT INTO mauer VALUES (1, ST_GeomFromText(
'POLYGON((-100 10, 100 10, 100 -10,
-100 -10, -100 10))', -1));
INSERT INTO segment VALUES (1, ST_GeomFromText(
'LINESTRING(-10 50, 10 -50)', -1));

Extracting the intersection between the rectangle and
the segment is simple and easy, as you already reported:

INSERT INTO segment_mauer
SELECT 1, ST_Intersection(mauer.geom, segment.geom)
FROM mauer, segment;

Final step: extracting the difference between the segment
and the polygon (or between the whole segment and the
shortest segment I've extracted just before, that is
the same in this specific case).

INSERT INTO segment_no_mauer
SELECT 1, ST_Difference(segment.geom, mauer.geom)
FROM mauer, segment;

INSERT INTO segment_no_mauer
SELECT 2, ST_Difference(segment.geom, segment_mauer.geom)
FROM segment_mauer, segment;

Both methods return exactly the same result, and it is
what we were expecting.
so, where is the problem ?

bye Sandro

[mj10777]

 Using your sample, there is no problem. The result are what was expected.

I have created an extract of the data for a 400x250 meters of the area shown in the first image.

There are 3 streets that cross over the Polygon.

When I adapt your script to also include the Polygon and the 3 streets

- the result is correct:

When called from my script, I get completely different results:

Yellow is the portion that should end at the top of the polygon.

Red is the correct result from the area inside the polygon

Green is the portion that should start at the bottom of the polygon.

2 Streets have 2 Linestrings, 1 only 1 Linestring (no green portion)

-- correct result:

SRID=3068;MULTILINESTRING((24810.59736691505 23664.22581131355,24819.85698453766 23644.2248388537),(24835.38784578578 23610.67784579614,24846.64252081577 23586.3675065906))

--incorrect result:
SRID=3068;MULTILINESTRING((24810.59736691505 23664.22581131355,24825.88165313041 23631.21142560983),(24825.88165313041 23631.21142560983,24846.64252081577 23586.36750659061))

Since the same command (ST_Difference) is being used in both cases, with the same data - but always with same results

- I fail to see what the cause is

Weird is also that the left line only returns a single LINESTRING, when the other 2 return MULTILINESTRINGs.

Is a puzzlement ...

The archive contains the extracted Database and the sql-scripts (yours - which has the added real Polygon/Linestring and mine).

It also contains a geotif of the map for that area and the QGIS project.

Mark

bye Sandro

[a.fu...@lqt.it]

Hi Mark,

Solved, mystery unveiled
caution: it may cause some severe headache !!!

short summary: any odd and puzzling result is simply caused by
the same numeric rounding issues that are plaguing the current
version of Topology.

recall: spatialite is based upon GEOS, and GEOS uses double
precision floating point coordinates that only supports
a finite precision. this can sometimes introduce nasty
truncation/rounding issues causing apparently crazy results.

you'll find a more thorough SQL analysis into the attached
SQL script, this is a quick summary:

1. the st_difference between soldner_segment and
soldner_segment_mauer always fails because the two
segments aren't really overlapping; a slight shift
affecting both start and end points causes an "X"
figure (by an incredibly small angle), but it's
enough to determine just a single intersection Point.

2. the st_difference between soldner_segment and
soldner_polygon correctly works.

3. a more sophisticated approach based on ST_Snap
(so to add more intersection Points) correctly
works in two cases but fails in a third case.

4. the final SQL query is useful in order to check the
actual presence of a real intersection between Start
End Points of soldner_segment_mauer and the second
geometry used by ST_Difference.

bye Sandro

[mj10777]

Thank you for the effort, I will look at the script tomorrow.

Yesterday I rewrote the script since the final task needed another approach that was not apparent in the first attempt

- I was expecting the same problem, but to my surprise it seems that all 306 result were done correctly

-- sometimes I have the impression that a step for step (instead of combinations of commands) somehow work differently

Mark

bye Sandro

[a.fu...@lqt.it]

a more robust approach:

1. precalculating all intersection points between polygon
rings and segments.
(temporarily saved as a single MultiPoint)

2. Snapping all segments against the above multipoint.
3. Snapping the polygon against the same multipoint.
these two steps together will ensure that now all
intersections are always explicitly represented by
some vertex, thus completely avoiding any furher need
to interpolate coordinates "on the fly" (the most
obvious cause introducing unexpected shifts).

4. as you can easily check, now ST_Difference correctly
works even when processing a couple of segments.

bye Sandro

[a.fu...@lqt.it]

On Thu, 12 Jan 2017 06:24:30 -0800 (PST), 'mj10777' via SpatiaLite
Users wrote:
> -- sometimes I have the impression that a step for step (instead of
> combinations of commands) somehow work differently
>

finite precision truncation/rounding issues happen in a rather
unpredictable way; they surely have a deteministic nature, but
can be easily affected by the exact order in which operations
are executed.
so it's not impossible that two apperently similar SQL queries
could be affected in differenct ways producing different results.

there is only one way safely ensuring a fair good consistency;
attempting to remove as far as possible any "flimsy" intersection
just based on points interpolated on the fly and very probably
affected by umpredictable slight shifts.

bye Sandro

[mj10777]

On Thursday, 12 January 2017 15:31:59 UTC+1, sandro furieri wrote:

a more robust approach:

1. precalculating all intersection points between polygon
    rings and segments.
    (temporarily saved as a single MultiPoint)

2. Snapping all segments against the above multipoint.

In my case I would 'export' the 306 rows into a work table, leaving the original alone 

3. Snapping the polygon against the same multipoint.
    these two steps together will ensure that now all
    intersections are always explicitly represented by
    some vertex, thus completely avoiding any furher need
    to interpolate coordinates "on the fly" (the most
    obvious cause introducing unexpected shifts).

Applying this on the original was something I was thinking about yesterday

- this method shows how it can be done