Search code examples
sql-serverpoint-in-polygonsqlgeography

Geography point within polygon showing incorrect results


In my database, I have a table containing a set of Polygons & their attributes and another table containing a set of points that have coordinates within the boundaries of those polygons. I've plugged QGIS into my database to ensure that they are located exactly where I expect them to be, and then ran a query to identify the polygon_id that each point is located within. This resulted in some very unexpected results, so I've done some further testing & tinkering and am still none the wiser.

I isolated one point and the polygon that it falls within, and this is how it's displayed by QGIS, using the same CRS (4326), connecting to my database:

enter image description here

However, when I use STContains or STIntersect I get a negative result.

When I run the following test the results suggest that the point does not exist with the polygon:

declare @p geography   
declare @l geography      

set @p = geography::STMPolyFromText('MULTIPOLYGON (((-2.01723932071019 57.1301507994778, -2.05959736945859 57.1308481812907, -2.06094305103518 57.1309872145762, -2.06442332315307 57.1311592624651, -2.06738702460041 57.131260961538, -2.06834761818132 57.1312257047487, -2.07352085461956 57.1308074773011, -2.08662715692557 57.1284264355531, -2.08835456332902 57.1282051063708, -2.0891720244275 57.1281861098907, -2.09070679418533 57.1283121889021, -2.09310621355588 57.1286955641499, -2.09502104689506 57.1293780321133, -2.09674004600622 57.1304564574772, -2.09753934746642 57.1297443781977, -2.09875197034386 57.1291864729083, -2.10340075289855 57.1282116756316, -2.10590378101433 57.1284283727703, -2.10916734028092 57.129042168269, -2.1099366825462 57.1291282326644, -2.11138826300974 57.1291886324809, -2.11294222137845 57.1287904597294, -2.11432776064015 57.1284126116958, -2.11590135588206 57.1272468512517, -2.11783267950974 57.1257414916239, -2.11972005111187 57.1230896037971, -2.12220448182387 57.1236735070728, -2.12237561434298 57.1239917795713, -2.12132109065263 57.1242488139312, -2.1173707565196 57.129607293332, -2.11658292484552 57.1328230202163, -2.11554089667428 57.1347555681128, -2.11504663847733 57.1359238965379, -2.11507076780594 57.1361851484318, -2.11527759405925 57.1364937439207, -2.11485673955606 57.1366369506911, -2.11452832733778 57.1372965941222, -2.1140394284489 57.1392757110769, -2.11377688986274 57.1399955310532, -2.11295601438074 57.1418667234446, -2.11208273346248 57.1427865720397, -2.10218925747153 57.1456872698938, -2.10042721504413 57.1441390669773, -2.09826717203099 57.1446680819023, -2.09796115562555 57.1448677424633, -2.09554608763705 57.1455415280263, -2.08978761299588 57.1409403132676, -2.08967555316753 57.1403717746151, -2.0882120280466 57.1405470942304, -2.07571111447462 57.1406560085666, -2.07235431383438 57.1418419845685, -2.06920884418592 57.1421474254686, -2.06521179825409 57.144067469068, -2.06100565229879 57.1456928396051, -2.01877021137073 57.1461574324006, -2.01746715794123 57.1443574286352, -2.01615980697691 57.1399056220016, -2.01632663387332 57.1358955754923, -2.01723932071019 57.1301507994778)))',4326);  
set @l = geography::STPointFromText('POINT (-2.0535344529468018 57.140210162062083)',4326);  

select @l.STContains(@p); 

I would expect the select statement to return a 1.

What's even stranger is that when I try and run that point against my entire polygon table:

select polygon_id from polygons where polygon.STcontains(@l)=1

It returns 246 results - out of 343 polygons.

I've tried both STContains, STWithin & STIntersects and all give similar results

I feel like I'm missing something very elementary, or QGIS is lying to me. Can anyone tell me where I'm going wrong?


Solution

  • SQL Server demands that polygons in a geography use Left-hand-Rule, which means that if you want to enclose a polygon you need to "walk" it with the inside of it on your left and the outside of it on your right.

    Your data is backwards, it's enclosing the outside of the polygon. You are also calling StContains on the wrong instance, you need to call it on @p and pass @l. You can reorient the polygon using ReorientObject()

    declare @p geography   
    declare @l geography      
    
    set @p = geography::STMPolyFromText('MULTIPOLYGON (((-2.01723932071019 57.1301507994778, -2.05959736945859 57.1308481812907, -2.06094305103518 57.1309872145762, -2.06442332315307 57.1311592624651, -2.06738702460041 57.131260961538, -2.06834761818132 57.1312257047487, -2.07352085461956 57.1308074773011, -2.08662715692557 57.1284264355531, -2.08835456332902 57.1282051063708, -2.0891720244275 57.1281861098907, -2.09070679418533 57.1283121889021, -2.09310621355588 57.1286955641499, -2.09502104689506 57.1293780321133, -2.09674004600622 57.1304564574772, -2.09753934746642 57.1297443781977, -2.09875197034386 57.1291864729083, -2.10340075289855 57.1282116756316, -2.10590378101433 57.1284283727703, -2.10916734028092 57.129042168269, -2.1099366825462 57.1291282326644, -2.11138826300974 57.1291886324809, -2.11294222137845 57.1287904597294, -2.11432776064015 57.1284126116958, -2.11590135588206 57.1272468512517, -2.11783267950974 57.1257414916239, -2.11972005111187 57.1230896037971, -2.12220448182387 57.1236735070728, -2.12237561434298 57.1239917795713, -2.12132109065263 57.1242488139312, -2.1173707565196 57.129607293332, -2.11658292484552 57.1328230202163, -2.11554089667428 57.1347555681128, -2.11504663847733 57.1359238965379, -2.11507076780594 57.1361851484318, -2.11527759405925 57.1364937439207, -2.11485673955606 57.1366369506911, -2.11452832733778 57.1372965941222, -2.1140394284489 57.1392757110769, -2.11377688986274 57.1399955310532, -2.11295601438074 57.1418667234446, -2.11208273346248 57.1427865720397, -2.10218925747153 57.1456872698938, -2.10042721504413 57.1441390669773, -2.09826717203099 57.1446680819023, -2.09796115562555 57.1448677424633, -2.09554608763705 57.1455415280263, -2.08978761299588 57.1409403132676, -2.08967555316753 57.1403717746151, -2.0882120280466 57.1405470942304, -2.07571111447462 57.1406560085666, -2.07235431383438 57.1418419845685, -2.06920884418592 57.1421474254686, -2.06521179825409 57.144067469068, -2.06100565229879 57.1456928396051, -2.01877021137073 57.1461574324006, -2.01746715794123 57.1443574286352, -2.01615980697691 57.1399056220016, -2.01632663387332 57.1358955754923, -2.01723932071019 57.1301507994778)))',4326);  
    set @l = geography::STPointFromText('POINT (-2.0535344529468018 57.140210162062083)',4326);  
    
    select @p.ReorientObject().STContains(@l);
    

    db<>fiddle

    To update your whole table you can do

    UPDATE YourTable
    SET ThePolygons = ThePolygons.ReorientObject()
    WHERE ThePolygons.STArea() >= 255030000000000;  -- more than half the world
    

    Note that if you actually have multiple polygons in the Multipolgon instances (which you don't here) then you'll need to break them out and work out if each one is oriented correctly.