SQL Server - Select all polygons nested within a larger polygon

Question

This is a seemingly simple SQL Server geometry question that I thought would have an out-of-the box solution, but I'm not having any luck finding one.

My intent is to select all the records within one table that have polygons that are nested (contained) within a larger polygon from another table. I had expected functions STWithin and STContains as being the solutions I needed, but unfortunately, both identify only the inner polygons within those nested in the larger polygon, not those nested polygons that are touching the boundary of the larger polygon. See image for example.

An alternative option that worked to my needs was STIntersection. The problem with this function, however, is that it only returns the geometry column! I would like to obtain the record ID instead. Does anyone have a suggestion on how this can be done?

STWithin:

select a.bg10 from
gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on a.shape.STWithin(b.shape) = 1
where b.mktname = 'Loop'

STContains:

select a.bg10 from
gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on b.shape.STContains(a.shape) = 1
where b.mktname = 'Loop'

STIntersection:

select a.shape.STIntersection(b.shape)
from gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on a.shape.STIntersects(b.shape) = 1
where b.mktname = 'Loop'

Edit:

One suggestion was to omit STIntersection and use solely STIntersects as follows:

STIntersects:

select a.bg10
from gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on a.shape.STIntersects(b.shape) = 1
where b.mktname = 'Loop'

The problem with this approach, is that STIntersects appears to select all polygons either within or outside and touching the larger polygon, not just those strictly within. See image for example.

Answer 1

In theory the queries that you have done should return the polygons you said haven't been returned. That makes me suspect that you might be encountering floating point error issues that SQL Server has with it spatial data types. Hence my comment about buffering the bounding polygon with a minimal amount.
So something like the following should get the results you want.

SELECT a.bg10 
FROM gis.usa_10_block_group a
    JOIN gis.usa_10_mkt_definition b
        ON a.shape.STWithin(b.shape.STBuffer(0.0001)) = 1
WHERE b.mktname = 'Loop'

Here's a quick example of the expected behaviour of a few of the spatial methods.

SELECT Geometry::STGeomFromText(WKT,0), Description
    , Geometry::STGeomFromText('POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))',0).STIntersects(Geometry::STGeomFromText(WKT,0)) Intersects
    , Geometry::STGeomFromText('POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))',0).STContains(Geometry::STGeomFromText(WKT,0)) Contained
    , Geometry::STGeomFromText('POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))',0).STOverlaps(Geometry::STGeomFromText(WKT,0)) Overlaps
    , Geometry::STGeomFromText('POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))',0).STTouches(Geometry::STGeomFromText(WKT,0)) Touches
FROM (VALUES
    ('POLYGON((0 0, 20 0, 20 20, 0 20, 0 0))'            ,'Interior corner')
    ,('POLYGON((90 90, 100 90, 100 100, 90 100, 90 90))' ,'Interior corner')
    ,('POLYGON((20 20, 40 20, 40 40, 20 40, 20 20))'     ,'Interior')
    ,('POLYGON((50 0, 70 0, 70 20, 50 20, 50 0))'        ,'Interior edge')
    ,('POLYGON((50 80, 70 80, 70 100, 50 100, 50 80))'   ,'Interior edge')
    ,('POLYGON((80 50, 100 50, 100 70, 80 70, 80 50))'   ,'Interior edge')
    ,('POLYGON((90 0, 110 0, 110 20, 90 20, 90 0))'      ,'Overlap')
    ,('POLYGON((100 50, 120 50, 120 70, 100 70, 100 50))','Exterior edge')
    )P(WKT,Description)
UNION ALL 
SELECT Geometry::STGeomFromText('POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))',0),'Bounding Area',null,null,null,null

Results

Description     Intersects Contained Overlaps Touches
--------------- ---------- --------- -------- -------
Interior corner 1          1         0        0
Interior corner 1          1         0        0
Interior        1          1         0        0
Interior edge   1          1         0        0
Interior edge   1          1         0        0
Interior edge   1          1         0        0
Overlap         1          0         1        0
Exterior edge   1          0         0        1
Bounding Area   NULL       NULL      NULL     NULL

Answer 2

The intersection query should look like this (assuming you want all records back from 'a'):

select a.* --get all columns from table 'a'
from gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on a.shape.STIntersects(b.shape) = 1
where b.mktname = 'Loop'

If you want just the areas of a that intersect b (ie. clipping a to b), you then add the STIntersection

select a.bg10
, a.STIntersection(b.geom) --clipped geometry from a against b
    from gis.usa_10_block_group a
    join gis.usa_10_mkt_definition b
    on a.shape.STIntersects(b.shape) = 1
    where b.mktname = 'Loop'

But this doesn't get you the polygons that are within b just yet...

This type of polygon-in-polygon is very testy with the boundaries and their conincidence - in order to be 'Within', the boundaries of a cannot be conincident with the boundaries of b - same goes for 'Contains'.

By these definitions, how many of your polygons in a are actually within b...?

So do you want to buffer b before you select polygons in a that are within? Or do a negative buffer on a?

Not sure what the exact answer here is...

Answer 3

As long as your scenario only has inner polygons that are either completely within or touching the boundary defined by the outer polygon then why not take the centroids of the inner polygons using EnvelopeCenter()? This is also much quicker because you are only asking whether a point is contained, rather than another polygon.

Answer 4

The query below works for two tables uploaded to mssql using ogr2ogr

SELECT * FROM table_with_features table_with_features INNER JOIN table_to_find_other_features_within table_to_find_other_features_within ON table_with_features.ogr_geometry.STIntersects(table_to_find_other_features_within.ogr_geometry) = 1 where table_to_find_other_features_within.attribute = 'whatever'