7

I am trying to select polygons that are within a larger polygon in PostGIS but those who share boundaries with the larger polygon are not being returned.

If I use ArcGIS Select By Location I get the result I need, 17 features. But if I use PostGIS I get only 4 features, the same as if I use ArcGIS "are completely within".

Correct result in ArcGIS

Queries that I have tried, but without success:

select count(a.cd_geocodi)
from t_ibge_setor_censitario a, t_ibge_municipio_2015 b
where b.nm_municip = 'PARAÍSO DAS ÁGUAS'

cand ST_Intersects(a.geom,b.geom); Result: 34

select count(a.cd_geocodi)
from t_ibge_setor_censitario a, t_ibge_municipio_2015 b 
where b.nm_municip = 'PARAÍSO DAS ÁGUAS'
and ST_Within(a.geom,b.geom);

Result: 4

select count(a.cd_geocodi)
from t_ibge_setor_censitario a, t_ibge_municipio_2015 b
where b.nm_municip = 'PARAÍSO DAS ÁGUAS'
and ST_Contains(b.geom,a.geom);

Result: 4

select count(a.cd_geocodi)
from t_ibge_setor_censitario a, t_ibge_municipio_2015 b
where b.nm_municip = 'PARAÍSO DAS ÁGUAS'
and ST_Coveredby(a.geom,b.geom);

Result: 4

select count(a.cd_geocodi)
from t_ibge_setor_censitario a, t_ibge_municipio_2015 b
where b.nm_municip = 'PARAÍSO DAS ÁGUAS'
and ST_Covers(b.geom,a.geom);

Result: 4

Is there a way to get all those 17 features in PostGIS?

  • whats your postgis query? try st_intersects, st_within,st_contains... – ziggy Feb 16 '17 at 15:28
  • ST_Intersects returns 34 features, because it returns features that share boundaries but is not within the features. ST_Within, ST_Contains, ST_Covers and ST_Coveredby always returns only 4 features. – Diego Henrique Feb 16 '17 at 15:41
  • 1
    ST_Within(ST_Centroid(b.geom), a.geom) – ziggy Feb 16 '17 at 16:00
  • 1
    You probably want ST_Intersects(a.geom, b.geom) AND NOT ST_Touches(a.geom, b.geom); I am assuming that the 34 results returned by the ST_Intersects are correct, but that you don't want ones that just share a boundary. – John Powell Feb 16 '17 at 16:50
  • 3
    It should be noted that ST_Within(ST_Centroid(b.geom), a.geom) is a horrible solution. It is possible for the centroid of a geometry to be outside the geometry (try a cresecent moon shape). It might work in this case, but ST_Intersects(a.geom, b.geom) AND NOT ST_Touches(a.geom, b.geom); is a much better solution. – John Powell Feb 16 '17 at 17:07
8

This depends a little on whether the small features are dependent on the larger features: do they actually share a common boundary, or are there small discrepancies along their apparently common edges?

ST_Within(a.geom, b.geom): gives you the geometries in b that are fully within a.geom. And it means fully within.

ST_Intersects(a.geom, b.geom): this is a little more tolerant, returning anything in b that merely touches a at all, even at a single point on the boundary. But if you combine this with NOT ST_Touches(a.geom, b.geom) as mentioned in a comment by John Barça (ST_Intersects(a.geom, b.geom) AND NOT ST_Touches(a.geom, b.geom)), then you will be left with (hopefully) the result you want. Note that small imperfections in the boundaries might still mean that the result isn't adequate, and you might need some kind of tolerance to get the correct result. If you still have issues, you could apply a small negative buffer to the geometries in b to make them smaller (this is pretty dirty though), or use ST_DWithin(a.geom, b.geom, 10), or you could perhaps ensure that a.geom and b.geom have a verified topological relationship and don't differ slightly at the edges.

  • This is a great answer. How would you verify a.geom and b.geom have a topological relationship? – ziggy Feb 17 '17 at 4:39
  • 3
    You might look at ST_CoveredBy, as well, to account for shared boundaries if topology is OK. It can avoid some subtleties of ST_Contains – Nate Wanner Feb 21 '17 at 18:21
0

I guess this is answer is a bit late, but I have figured out a way to do it in QGIS, and a similar process can perhaps be used for PostGIS. I first used the geoprocessing, buffer tool. Set very small values (this will depend on the scale of your features) to get a shape file having very similar boundaries.

After that, I used the new layer as the basis for selecting by location, using within.

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.