1

I have a point layer (with a buffer value), and a polygon layer (footprint) that intersects the point buffer.

Using only SQL in QGIS Virtual Layers (i.e. based on Spatialite), how can I get the difference between the point buffer (buffer(pt.geometry,bufvalue)) and footprint, where the difference does not intersect with the point geometry?

See the following image for an illustration - green is the footprint, orange pentagon is the point, orange polygon is the difference. I just want the ticked geometry.

enter image description here

The difference returns a multipolygon, and I only need one part. However, I can't split multipolygons into single polygons (there is no ST_Dump() equivalent in Spatialite due to SQLite not handling array values, apparently)

And I can't use the CLI/GUI commands (which separate multigeometries into 'elementary geometries') as I want to create a virtual layer in QGIS, so it has to be 100% in SQL.

How can I approach this problem using SQL in Virtual Layers only? (i.e. not using any other processing tools)

I am thinking I may have to loop over the number of geometries for multipolygons using recursive functions in place of generate_series(), and test if each one intersects with the point geometry...

  • SpatiaLite queries can be combined with QGIS expressions when creating the virtual layer. For example, to keep the indexed part as 1, the following clause can be used: SELECT geometry_n(geometry,1) as geom FROM 'layer'. Therefore, if you could know the index of the part to be extracted (even the parts should be able to be sorted with order_parts()) I think the problem would be reduced to how to implement the expressions in the query. – Gabriel De Luca Nov 9 '19 at 8:34
  • Unfortunately I have no way of directly knowing what the index of the part is to be extracted - each part has to be spatially tested (disjoint from point geometry). The issue is that I can't test each part separately without 'splitting' the multipolygon into single parts, which I can't do with Spatialite SQL. So I have resorted to my convoluted fallback of looping through each part (where num_geometries(geometry) >1) using a sequence of indices generated by a CTE in geometry_n(), then keeping the ones that don't intersect with the point geometry... then collecting those back... cripes. – she_weeds Nov 9 '19 at 10:33
  • 1
    Today I also wrote a query to dump geometries in a virtual layer, for the following question: gis.stackexchange.com/questions/334371/… . If you wanted to see it and give your opinion, the same if you wanted to write in an answer how you achieved your index sequence, I think it would be a good contribution. – Gabriel De Luca Nov 12 '19 at 16:58
1

I was trying to improve my query a bit, to dump single-part geometries in a virtual layer. I would appreciate if you can test it with your data:

WITH RECURSIVE parts AS (
  SELECT d.buffer_fid, d.poly_fid, d.total, 1 AS part
  FROM dump_me AS d
  UNION ALL
    SELECT d.buffer_fid, d.poly_fid, d.total, part + 1 AS part
    FROM dump_me AS d 
    INNER JOIN parts as p
      ON d.buffer_fid = p.buffer_fid AND d.poly_fid = p.poly_fid AND p.part < d.total
),
dump_me AS (
  SELECT buffer_fid, poly_fid, NumGeometries( geom) AS total
  FROM diff
),
diff AS (
  SELECT b.fid AS buffer_fid, p.fid AS poly_fid , ST_Difference( b.geom, p.geometry) AS geom
  FROM buffers AS b
  INNER JOIN polys AS p
    ON ST_Intersects( b.geom, p.geometry) = 1
),
buffers AS (
  SELECT fid, ST_Buffer( geometry, <bufvalue>) AS geom
  FROM points
),
dumped AS (
SELECT p.buffer_fid, p.poly_fid, p.part, GeometryN( d.geom, p.part) as geom
FROM parts AS p
INNER JOIN diff AS d
  ON p.buffer_fid = d.buffer_fid AND p.poly_fid = d.poly_fid
)
--SELECT * FROM buffers;
--SELECT * FROM diff;
--SELECT * FROM dump_me;
--SELECT * FROM parts;
--SELECT * FROM dumped;
-- to debug, comment on the following lines and uncomment the previous SELECT clause of your interest
SELECT d.buffer_fid, d.poly_fid, d.part, d.geom
FROM dumped AS d
INNER JOIN points AS p
 ON  d.buffer_fid = p.fid AND ST_Intersects( d.geom, p.geometry) = 0;

I am assuming the existence of a points layer and the existence of a polys layer, both with unique values for the fid field.

The same would appreciate the opinion of any SQL specialist to improve the query using tools available in SpatiaLite, in the virtual layers environment of QGIS.

  • Thanks Gabriel, excellent answer! This by and large worked - it is quite similar to the query I used except much cleaner :-) and I kept my recursive CTE to a hard limit of 5 in order to minimise processing (which works in my use case knowing the data). For anyone else bumbling their way through I realised after much trial and error that the recursive CTE must be at the top for it to work. – she_weeds Nov 25 '19 at 8:39
  • @she_weeds You are welcome. It was interesting for me to also understand how the recursive CTE works. – Gabriel De Luca Nov 25 '19 at 11:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.