# Determining if polygon not enclosed by other polygons [closed]

Does anyone know of a routine/function or transformer that would identify if a polygon is not fully enclosed/surrounded by other polygons?

The polygons have no gaps or slivers between them so if a certain section of polygon is not attached to another polygon then these are the ones I want to keep.

I am happy to use any of PostGIS/FME/QGIS I just cannot seem to find a function or transformer that would do it.

I have added a picture to help explain:

I am after the purple polygons and ignore the yellow ones as they are not fully enclosed by other squares

• Are the polygons on the same layer? And if you worked with USA data then Texas would be selected but Colorado not, or? Commented Aug 6, 2018 at 9:59
• Yes the same layer. Thats exactly what I am after. Commented Aug 6, 2018 at 10:44
• I undestood that if a certain section of polygon is not attached to another polygon then these are the ones I want to keep means the purple polygons. Commented Aug 6, 2018 at 11:05
• sorry - yes your are right, i just edited. Mondays! Commented Aug 6, 2018 at 11:11
• So your question is "How can I tell if a feature in a set of features has any external boundaries?". Test the linear intersection of A with the dissolved features. Commented Aug 6, 2018 at 13:18

Assuming your polygons are topologically correct (no gaps, or overlaps), then a polygon A is enclosed by other polygons if its perimeter is equal to the sum of intersection length between its boundary and the boundary of other polygons. This method lends itself to use of a spatial index and avoids performing a relatively costly union of all other polygons.

SELECT a.id
FROM my_data a
INNER JOIN my_data b ON (ST_Intersects(a.geom, b.geom) AND a.id != b.id)
GROUP BY a.id
HAVING 1e-6 >
abs(ST_Length(ST_ExteriorRing(a.geom)) -
sum(ST_Length(ST_Intersection(ST_Exteriorring(a.geom), ST_ExteriorRing(b.geom)))));


I have an answer that should be logically correct but the SQL part is not complete. The idea is based on intersections of adjacent polygons. The intesections are either lines if polygons share segment(s), or points if polygons meet at one vertex.

The first step is to compute the intersection of one polygon and the union of all the polygons on the layer, except the selected polygon.

The intersection of polygon 1 and union of all the rest is a multilinestring with two members. The intersection of polygon 6 and anything else but polygon 6 consists of four linestrings members which are connected.

The second step is to try if it is possible to build a new polygon from the intersections. For the connected lines from polygon 6 it is possible but not for the two linestrings from polygon 1.

This is the SQL that I used for testing.

select ST_MakePolygon(ST_LineMerge(c.intersection)) from
(select
ST_Intersection(sub_a.a,sub_b.b) as intersection from
(select ST_Union("GEOMETRY") as a from outline_test where gid!=6) as sub_a,
(select "GEOMETRY" as b from outline_test where gid=6) as sub_b) as c;


With gid=1 the result is an error:

ERROR:  lwpoly_from_lwlines: shell must have at least 4 points


and with gid=2:

ERROR:  lwpoly_from_lwlines: shell must be closed


Another approach: The polygons you want to find have 1-dimensional intersection (lines) with the outer ring of the whole area as an union. If polygon is only touching the outer boundary the intersection is a point with dimension=0, and polygons which are totally inside to not intersect as all (so they are disjoint).

So you can select the features with ST_Dimension=1 from SQL query like

select outline_test.gid,ST_Dimension
(ST_Intersection(outline_test."GEOMETRY",sub_a.geometry))
from
(select ST_ExteriorRing(ST_Union("GEOMETRY")) as geometry from outline_test)
as sub_a, outline_test;


The query does not make difference between "touches" and "disjoint" as documented in https://postgis.net/docs/ST_Dimension.html If the dimension is unknown (empty GEOMETRYCOLLECTION) 0 is returned but that does not matter in your use case. Rereffing to the data and image in my other answer, other polygons than numbers 5, 6, 8, and 10 should be found and that happens.

1;1
2;1
3;1
4;1
5;0
6;0
7;1
8;0
9;1
10;0
11;1
12;1
13;1
14;1


I am sure it is possible to write more elegant queries based on the outer boundary and DE-9IM relations http://postgis.net/docs/using_postgis_dbmanagement.html#DE-9IM.

• I think you've expressed the inverse problem with "The polygons you want to find have 1-dimensional intersection (lines) with the outer ring of the whole area as an union" - the required polygons are the yellow ones in the diagram and they don't have 1-d intersection with the external. Commented Aug 6, 2018 at 15:45
• OP wants the purple polygons. Commented Aug 6, 2018 at 15:57
• Oh then the title is misleading: "Determine if a polygon is enclosed by other polygons"... I'll edit.. Commented Aug 6, 2018 at 17:04

There is a solution in R, using spdep::poly2nb. This tools checks for the two types of neighborhood: queen and rook, and also includes snap feature to account for topographic inconsistencies/slivers or gaps.

Here is a reproducible example:

library(spdep)    # neighbours
library(raster)

r <- raster(nrow=11, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(NA,  NA, NA, NA, NA,1,
NA, NA, NA, 1, 1, 1,
NA, NA, 2, 1, 3, 1,
NA, NA, 1, 1, 1, 1,
NA, NA, 1, 2, 2, 1,
NA, NA, 1, 1, 1, NA,
NA, 1, 1, 1, 1, NA,
NA, 1, 1, 1, 1, NA,
NA, 1, 1, 1, 1, NA,
NA, 1, 1, 1, NA, NA,
NA, 1, NA, NA, NA, NA),
nrow = 11,
ncol = 6,
byrow = TRUE)

# Convert raster to polygon
# creates SpatialPolygonsDataFrame
polys <- rasterToPolygons(r)


Identify if the polygons are on the edge based on nearest neighbors. The neighborhood can change based on queen = TRUE or queen = FALSE, which equals to rook neighborhood.

nb <- poly2nb(polys,
queen=TRUE, # queen = FALSE refferes to Rook neighborhood
#row.names = polys,
snap = 0) # snap corrects for the gaps/slivers


Get the number of neighbors using card() function, and include this as and attribute to SpatialPolygonDataFrame:

polys$nb_count<- card(nb) # Has the stand an open edge? Is surrounded by neighbors, pre-value is FALSE polys$open_edge = ifelse(card(nb) <max(card(nb)), 1, 0)

# If has complete neighbors,
# check the differences in height
windows()
par(mfrow=c(2,1))
spplot(polys,
"open_edge") # plot value 'open_edge' 1 or 0


Which result:

It's easy to get it,

if I have understood the translation correctly by looking at the picture you presented.

Input data is shown in the figure

Run the script:

WITH tbla AS (SELECT ST_Union(geom) geom FROM grid),
tblb AS (SELECT ST_ExteriorRing(ST_Buffer((ST_Buffer((geom),-0.0001)), 0.0001)) geom FROM tbla)
SELECT (a.geom) geom FROM grid a JOIN tblb b ON ST_Intersects(a.geom,b.geom)


Look at the result

Tested by: PostgreSQL 11.1, PostGIS 2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1...

For your example, create a geo function:

CREATE OR REPLACE FUNCTION ST_ExtractionExternalGridCellsByDoubleBufferization(
geom GEOMETRY,
RETURNS TABLE (geom GEOMETRY) AS
$BODY$
WITH
tbla AS (SELECT (ST_Dump($1)).geom), tblb AS (SELECT ST_ExteriorRing(ST_Buffer((ST_Buffer(ST_Union(geom),-radius_buffer)), radius_buffer)) geom FROM tbla) SELECT ST_Collect(a.geom) geom FROM tbla a JOIN tblb b ON ST_Intersects(a.geom, b.geom);$BODY\$
LANGUAGE SQL


RUN

SELECT ST_ExtractionExternalGridCellsByDoubleBufferization(ST_Collect(geom), 0.001) geom FROM grid


So, in FME you would probably use the NeighborFinder and count the number of neighbors. If you included diagonal polygons, then you just need to check if there are 8 neighbors, but since you don't seem to need that, we need to check for 4 neighbors, but only horizontally/vertically.

So, here's a workspace to do that:

Basically:

• Use a CentrePointReplacer to convert polygons to centre points
• Use a NeighborFinder:
• Set Candidate Only Mode
• Neighbors to Find = 4
• Maximum distance = width/height of your polygons
• Closest Candidate List Name = list
• Use a ListElementCounter to count the number of neighbors found
• Use a Tester to test if _element_count < 4
• Use a SpatialFilter to filter out the original data
• Tester:Passed = Filter, Original Data = Candidates
• Predicate to Test = "Filter is Within Candidate"

The SpatialFilter:Passed output is the data you need.

I put my demo workspace on Dropbox if you want to check it out.