ogr2ogr sql - returning points from one sqlite dataset that intersect with another

Question

I am having great difficulty with the simple st_intersect command when applied to two datasets containing points, derived from the same source. I am using osgeo4w gdal 2.1.3 at command line.

My aim is simply to return in the output the point geometry features from sqlite table A where they intersect with sqlite table B (also point geometry). There are no other geometries present in the datasets. Whilst the command seems to work ok for most of the data (providing the expected output), for a small percentage of the data the st_intersect will not identify points which clearly do intersect upon visual scrutiny - i.e. after zooming to 1:1! These points do not appear in the output. I am puzzled. Might there be a precision issue here re. point geometry location, or is there some other reason for the omission of points? What method will achieve the desired result?

ogr2ogr -unsetFid -nlt POINT -f SQLite -dsco SPATIALITE=YES outfile.sqlite infile.sqlite -nln outtable -dialect sqlite -sql "SELECT A.* FROM 'intable1' A, 'intable2' B WHERE ST_Intersects(A.geometry, B.geometry)"

Answer 1

I solved this with a union/ snap to grid query, performed on each input table in the initial question. The intersect then behaved as expected on the output tables although output point locations where very slightly misaligned on occasion. Though not ideal from the point of view of visual analysis, this met my requirements as a solution.

ogr2ogr -unsetFid -nlt MULTIPOINT -f SQLite -dsco SPATIALITE=YES \
   outfilesnap.sqlite infile.sqlite -nln outtablesnap1 -dialect sqlite \
   -sql "SELECT ST_Union(ST_SnapToGrid(geometry,0.00000006666)),ogc_fid \
   FROM 'intablesnap1' GROUP BY ogc_fid"

ogr2ogr -unsetFid -nlt MULTIPOINT -f SQLite -dsco SPATIALITE=YES \
   outfilesnap.sqlite infile.sqlite -nln outtablesnap2 -dialect sqlite \ 
   -sql "SELECT ST_Union(ST_SnapToGrid(geometry,0.00000006666)),ogc_fid \
   FROM 'intablesnap2' GROUP BY ogc_fid"

ogr2ogr -unsetFid -nlt MULTIPOINT -f SQLite -dsco SPATIALITE=YES \
   outfileintersect.sqlite outfilesnap.sqlite -nln outtableintersect \
   -dialect sqlite -sql "SELECT A.* FROM 'outtablesnap1' A, 'outtablesnap2' B \ 
   WHERE ST_Intersects(A.geometry, B.geometry)"

I also found that altering the intersect sql query to "SELECT * FROM 'outtablesnap1' WHERE geometry NOT IN (SELECT A.geometry FROM 'outtablesnap1' A, 'outtablesnap2' B WHERE ST_Intersects(A.geometry, B.geometry))" returned the inverse (i.e. points in the first table that did not intersect with the second table). This was most useful.