polygon intersection and null geometry

[romain riviere]

Hi,

I'm a new spatialite user and I'm french.

I've tried , trought spatialite-gui, to create an intersection table from two polygons tables with that query:

create table "table3" as select

st_multi( st_intersection( table1.geometry,table2.geometry) ) as geometry

from table1, table2

where

st_intersects( table1.geometry,table2.geometry )

But when I do recovergeometrycolumn() with MULTIPOLYGON, the query failed (constraint failed). I think this is because the geometry column contains NULL geometry.

i don't understand why the query leads to NULL geometry: st_multi and intersects() should avoid it, no ?

Best regards

Romain,

[a.furieri]

Hi Romain,

I've performed a simple test using
the ISTAT-2010 Counties and Regions
(the same used in the Cookbook,
recipe #19)

CREATE TABLE test AS
SELECT c.nome_pro AS nome_pro, r.nome_reg AS nome_reg,
ST_Multi( ST_Intersection( c.geometry, r.geometry) ) as geometry
FROM counties AS c, regions AS r
WHERE ST_Intersects( c.geometry, r.geometry );

SELECT DISTINCT ST_GeometryType(geometry),
ST_Srid(geometry) FROM test;

> GEOMETRYCOLLECTION 23032
> MULTILINESTRING 23032
> MULTIPOLYGON 23032

SELECT nome_pro, nome_reg, geometry,
ST_GeometryType(geometry)
FROM test
WHERE ST_GeometryType(geometry) =
'MULTILINESTRING' OR
ST_GeometryType(geometry) =
'GEOMETRYCOLLECTION';

===

Explanation:
a) any MULTILINESTRING represents a
boundary shared by a County and
the adjacent Region
b) GEOMETRYCOLLECTION is more or
less the same, but this time
there is some sparse POINT as well.

Conclusion:
ST_Intersects() computes *any*
intersection, this including
segments and points.

So rewriting the query as follows
seems to be a better solution:

CREATE TABLE test2 AS
SELECT c.nome_pro AS nome_pro, r.nome_reg AS nome_reg,
ST_Multi( ST_Intersection( c.geometry, r.geometry) ) as geometry
FROM counties AS c, regions AS r
WHERE ST_Intersects( c.geometry, r.geometry )
AND ST_GeometryType ( ST_Multi(
ST_Intersection( c.geometry, r.geometry) )) =
'MULTIPOLYGON';

Alternatively, we can use the first query, and
then delete any 'odd' result:

DELETE FROM test
WHERE ST_GeometryType(geometry) <>
'MULTIPOLYGON';

bye
Sandro

[markb]

Romain,
You are correct that the intersection of two polygons should only result in another polygon (or multipolygon).  And your ST_Multi CAST function should solve the single/multi constraint validation, but not as you can tell with a NULL geometry in the mix.  My only idea, assuming problem free intersects/intersection functioning, is that you might have a NULL geometry in one of your source polygon tables (maybe the result of a prior overlay operation in your workflow, particularly on an un-recovered geom field, but not limited to it).  The f(x) "Intersects" can return a NULL value, above and beyond just Boolean, though this is a -1 (not really a NULL).  SQLite rejects only zero/false/actual NULL values in any of its Boolean operations (e.g. WHERE clauses), so all other values pass as true and remain unfiltered.  If you think it is a NULL problem, try to pinpoint it.  Here are some pertinent NULL links: 

SpatiaLite 2.4.0rc4 Intersects f(x) -
http://www.gaia-gis.it/spatialite-2.4.0-4/spatialite-sql-2.4-4.html#p12

SQLite Boolean Expressions -
http://www.sqlite.org/lang_expr.html#booleanexpr

Hope this helps, good luck.  -Mark

[markb]

I stand corrected, I can't believe I've never encountered this scenario.  I somehow thought dimensionality was maintained.  Go figure, you live you learn.  Obviously still breaking-down ESRI centric geoprocessing thinking, even after all these years.  So another variant, to play off of Sandro's example, this one only uses spatial predicates (Boolean returns) in the WHERE clause: 

CREATE TABLE test2 AS
SELECT c.nome_pro AS nome_pro, r.nome_reg AS nome_reg,
 ST_Multi( ST_Intersection( c.geometry, r.geometry) ) as geometry
FROM counties AS c, regions AS r

WHERE ST_Overlaps( c.geometry, r.geometry )
  OR ST_Equals( c.geometry, r.geometry);

[romain]

Thank you fro your help.

A geometry check of my two source tables only shows "polygon", not a
single "NULL" value.
In the where clause, I've tried (overlaps or equals) but the result
still contains a null value.
Nevertheless, adding the st_multi()="multipolygon" works fine !

I've tried somethink very strange (to my mind):
1/ Create table3 (.....)
2/ addgeometrycolumn()
3/ insert into table3
select
st_multi(st_intersection(table1.geometry,table2.geometry) as geometry
where
ST_Intersects(table1.geometry,table2.geometry)

I still have the null value, but the geometry column's update works
fine !

So I think there are 3 solutions:
1/ create first a table, a geometry column and then insert datas
2/ Complete the SQL clause with
st_geometrytype(st_multi(st_intersection(...))="multipolygon"
-> But, this condition may slow down the query execution, am I
right ?
3/ Delete the 'odd' result by hand before recovering geometry column

[markb]

I'm afraid to use the word "definitely" anymore, but it seems something is off.  Have you had a chance to inspect your NULL output record, specifically going back to its parents?  If you don't have a malformed geometry (checked with IsValid), then something else appears to be occurring. 

This "could" be an alternative fix (or diagnosis tool), without the geom create penalty: 

CREATE TABLE test2 AS
SELECT c.nome_pro AS nome_pro, r.nome_reg AS nome_reg,

 ST_Multi(ST_Intersection(c.geometry, r.geometry)) as geometry

FROM counties AS c, regions AS r
WHERE

 CASE ST_Overlaps(c.geometry, r.geometry)
  WHEN 0 THEN 0
  WHEN 1 THEN 1
  ELSE 0 --changes -1/unknown to 0/false
 END
 OR
 CASE ST_Equals(c.geometry, r.geometry)
  WHEN 0 THEN 0
  WHEN 1 THEN 1
  ELSE 0 --changes -1/unknown to 0/false
 END;

[romain]

I've check geometries from my two source tables, (gui->malformed
geometry)
and it appears that the first one ( Regions from test-2.3.sqlite)
contains many invalid geometries:

I've 14 "repeated vertex" marked "reapairable" and 1 "GEOS invalid
geometry".

This seems to be the reason of my issue....

Thank you so much !

So, it seems that, in order to use spatialite correctly, valid data is
required.... but with mapinfo/ESRI users...this type of data doesn't
exist !

[markb]

Glad to help, your mapinfo/ESRI quip was hilarious.  

[markb]

I see some potential pitfalls with the assistance we've given you; I'm still adjusting my worldview regarding this discussion.  It is more than conceivable, now that I'm fully considering the complexities, that even though the WHERE clauses we've supplied you guarantee poly(s) for each output record, that each will have certain significant shortcomings.  And this also uncovers another more problematic issue...

My method could result in several geometry types, since Overlaps OR Equals usage doesn't exclude points or lines; it just makes sure polygons are present, but not exclusively.   You could end up with a GeometryCollection being returned by your Intersection (I've tested & verified this).  Now for the worse news.  ST_Multi doesn’t remove points or lines, just Multi wraps them.  So a mixed type GeomColl would stay a GeomColl.  CastToMultiPolygon doesn't work on GeometryCollections with mixed sub-geometries (not sure about Polygons with MultiPolygons?) as per Sandro's documentation (I double checked this), since it returns NULL in those instances (i.e. Poly with Pt and/or Line).  No methods exist to extract specific geom-types from mixed GeomColls, so your output remains trapped in those instances and can't be separated and pooled with others of the same type.  Related to the very reason Sandro created the geom-type casting functions to begin with.  Though I haven't encountered this myself (beyond testing here & now), I can see how some datasets would be more vulnerable to this situation than others (I've just been lucky).  This eventuality is quite unsettling, though understandingly easy to miss. 

Sandro's example definitely yields exclusively MultiPolygons, but would eliminate any polygons contained within a GeomColl (by result of points and/or lines also).  This could result in significant omissions.  Also very easy to overlook, particularly with topologicaly nice political boundaries (verses natural phenomena).  This leads back to retaining GeomColls and the trappings currently entailed by that, as stated above. 

Intersection works the way it does for a reason; one that I understand now (and agree with), but it definitely requires certain considerations (particularly with multiple intersections per single part geoms and even worse with multiparts on multiparts).  One has to plan their path accordingly and will need the right tools to navigate through it.  Maybe ExtractTo* or modify CastTo* (though possibly misleading).  I hope my testing is more of a help than a bother, but either way sorry for another feature request Sandro. 

Thanks for your patience as always everybody.  -Mark

[a.furieri]

Hi Mark,

all this make perfect sense to me.

SELECT ST_AsText(ExtractMultiPoint(ST_GeomFromText(
'GEOMETRYCOLLECTION(POINT(1 1), LINESTRING(4 5, 6 7, 8 9), POINT(30
30))' )));
> MULTIPOINT(1 1, 30 30)

SELECT ST_AsText(ExtractMultiLinestring(ST_GeomFromText(
'GEOMETRYCOLLECTION(POINT(1 1), LINESTRING(4 5, 6 7, 8 9), POINT(30
30))' )));
> MULTILINESTRING((4 5, 6 7, 8 9))

SELECT ST_AsText(ExtractMultiPolygon(ST_GeomFromText(
'GEOMETRYCOLLECTION(POINT(1 1), LINESTRING(4 5, 6 7, 8 9), POINT(30
30))' )));
> NULL

SELECT ST_AsText(ExtractMultiPoint(ST_GeomFromText(
'MULTIPOINT(1 1, 30 30)' )));
> MULTIPOINT(1 1, 30 30)

SELECT ST_AsText(ExtractMultiLinestring(ST_GeomFromText(
'LINESTRING(4 5, 6 7, 8 9)' )));
> MULTILINESTRING((4 5, 6 7, 8 9))

... and so on ...

===================
already committed to SVN [rev.106]:
svn https://www.gaia-gis.it/svn/libspatialite libspatialite-svn

bye Sandro

[markb]

Eccellente, thank you so much Sandro.  Now I can separate my spaghetti sauce (polys) from my linguine (lines) and meatballs (pts) diner, a la SpatiaLite Cookbook style.  So my final answer to obtain all areas-in-common (i.e. polygonal only intersection) between two polygon tables is as follows: 

CREATE TABLE test2 AS
SELECT c.nome_pro AS nome_pro, r.nome_reg AS nome_reg,

ExtractMultiPolygon(ST_Intersection(c.geometry,r.geometry)) AS geometry

FROM counties AS c, regions AS r

WHERE ST_Overlaps(c.geometry, r.geometry)
  OR ST_Equals(c.geometry, r.geometry);

This is functionally about as concise as it gets I think, and should be equivalent to an ESRI style intersection (for those migrating from ESRI to the wonderful SpatiaLite).  -Mark