3

I get different results, performing a spatial join with PostGIS and QGIS.
The data can be downloaded here.

I try to join vg2500_krs (402 rows) with vg2500_bld (16 rows). I am interested in writing the information from vg2500_bld into vg2500_krs, keeping all records from vg2500_krs.
In QGIS I select Vector - Data Management Tools - Join Attributes by Location. The resulting shapefile has 402 rows, as expected.

However, when I try to do this with PostGIS, I get a shapefile with 638 rows.

This is the query I use:

CREATE TABLE krs_bld AS
SELECT k.gen, k.rs, k.geom, b.name_bld
FROM vg2500_krs k
LEFT OUTER JOIN vg2500_bld b
ON ST_Intersects(k.geom, b.geom) AND NOT ST_Touches(k.geom, b.geom);

Further, when I look at the Shapefile created by PostGIS in QGIS, especially at the duplicated rows, they "do not exists". When I highlight them and "Zoom to selected Feature", there is nothing highlighted in the mapview. It is as if nothing was selected. But, when I calculate the area for the rows, it does not result in 0. Where does this behavior come from?

2
  • For the on statement you could try: ON ST_Intersects(ST_Centroid(k.geom), b.geom) Nov 30, 2016 at 11:17
  • @Vesanto that does not explain, why it does not "respect" the ST_Touches...
    – four-eyes
    Nov 30, 2016 at 11:28

1 Answer 1

2

You can remove the ST_Touches without impacting the results of the SQL query.

AND NOT ST_Touches(k.geom, b.geom);

The reason for this that ST_Touches and PostGIS in general is really precise. To you it looks like the two datasets line up perfectly. But they do not.

If you really zoom in on the two files at the boundaries:

enter image description here

They do not line up. ST_Touches works on a perfect touch, this is a precision you cannot achieve in QGIS or ArcGIS, you would need a true topological GIS (like GRASS).

So the workaround, and how QGIS is probably already doing it, is to use the centroids for the intersection. Or even better a point on surface.

So if you modify the query to:

CREATE TABLE krs_bld AS
SELECT k.gen, k.rs, k.geom, b.name_bld
FROM vg2500_krs k
LEFT OUTER JOIN vg2500_bld b
ON ST_Intersects(ST_PointOnSurface(k.geom), b.geom);

You should get the results you desire.

You can read about a similar issue here, and why ST_Centroid might not work either.

1
  • Oh wow. That is something I did not know! Thanks a lot for pointing that one out. You taught me something
    – four-eyes
    Dec 1, 2016 at 17:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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