I have two tables, 'protein' and 'data_aggreg', in a PostGIS database. I am simply trying to query the points that intersect (all of the points in 'data_aggreg' share exactly the same geometry as a portion of the points in 'protein' table).

This is typically a simple operation. However, this time PostGIS and QGIS return nothing and ArcGIS returns the correct selection of points. The data in both QGIS and ArcGIS are sourced directly from the same PostGIS database via their respective methods of connecting to external DBs (ie, no shapefiles, csv's, etc are being used for this analysis). Both tables are in projected coordinate system EPSG 26912 and both tables are single point geometry.

When I run ST_Intersects in PostGIS on the two geom columns nothing is returned. When I run the 'select by location' tool in QGIS nothing is returned (I've also tried spatial join with same effect).

Syntax for ST_Intersects, which returns nothing:

SELECT protein.protein, data_aggreg.pro_2016
from production.protein, scratch.data_aggreg
where ST_Intersects(geom26912, geom);

However when I run 'select by location' in ArcGIS I get the expected result. That is, the subset of points in the 'protein' table that share the same geometry with 'data_aggreg' are correctly selected. Also, when I change my PostGIS 'where' query from ST_Intersects to geom1 = geom2 I also get the expected result.

Syntax when using '=' in PostGIS, which returns the expected result:

SELECT protein.protein, data_aggreg.pro_2016
from production.protein, scratch.data_aggreg
where geom26912 = geom;

It is my understanding that if any part of features share any space then ST_Intersects will return 'true'. So it follows that if geom1 = geom2 then ST_Intersects should return 'true'. But this is not the case.

So why the difference in results within PostGIS when using ST_Intersects vs '=', and why the difference among PostGIS when using ST_Intersects, QGIS 'select by location', and ArcGIS 'select by location'?

  • It may be helpful to post the syntax you used for postgres.
    – Aaron
    Nov 16 '16 at 18:15
  • Both SELECT ST_Intersects('POINT(0 0)'::geometry, 'POINT(0 0)'::geometry); and SELECT 'POINT(0 0)'::geometry = 'POINT(0 0)'::geometry; return true in my tests ...
    – underdark
    Nov 16 '16 at 18:18
  • Shouldn't you still use Equals when comparing geometries instead of "="?
    – user30184
    Nov 16 '16 at 18:34
  • What if you add ST_SnapToGrid with a rather small tolerance to ST_Intersects?
    – user30184
    Nov 16 '16 at 18:37
  • 1
    I have a feeling that the problem is in the main SQL but I must think a bit further.
    – user30184
    Nov 16 '16 at 18:47

= in PostGIS means bounding box equality. There is long history of that and people have debated changing it. But that's where that stands.

So for anything but points (even points you can't quite trust because of box rounding), it does not mean what you think it means.

As someone mentioned use something like ST_Equals (note this doesn't care about the orientation geometries), or ST_OrderingEquals.

I think ST_OrderingEquals is more akin to what people think of when they think of equality in the database sense.

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.