2

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)"

  • 1
    Intersect for points means the same as equals and coordinates must be the same till the last digit. Perhaps introducing some tolerance with SnapToGrid helps. Or build a buffer for either a.geometry or b.geometry. – user30184 Apr 24 '17 at 19:39
0

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.

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.