STIntersect says every point is intersecting every zone - need to confirm

Question

I am implementing a Vehicle Safety program for an industry group. The group have provided 3k+ geozones that cover sites that are members of the group. The file provided was a shapefile. My colleague extracted it and converted it to 4326 and uploaded it into a new database. I checked the zones using IsValidDetailed(). 7 were invalid, but righted themselves when MakeValid() was run on them. The geography table SaferTogether has had a spatial index created, just using auto.

I was initially testing using the following. This would run every point against every zone, and I got results after just over an hour.

INSERT INTO @Intersects (LogId, ogr_fid, NodeId)
SELECT p.LogId, s.ogr_fid, p.NodeId
FROM LogPos p,
     [gis]..safertogether as s
WHERE p.MyPoint.STIntersects(s.geog4326) = 1
and AssembledTime between @PeriodStart and @PeriodEnd

I'm now testing a single trip of 27 log points, and I've noticed I have 81162 results This roughly matches the 27 points x 3k+ zones. I don't believe that every zone overlays every point. The code below is the testing for one trip worth of logs. It is not efficient.

UPDATE #Logs
SET IsIntersect = 0

SET @Message = 'Detecting Log intersects with Zone ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogInfo, @Section

UPDATE #Logs
SET IsIntersect = 1
WHERE @Geography.STIntersects(MyPoint) = 1

SELECT @rc = COUNT(1)
FROM #Logs WHERE IsIntersect = 1

SET @Message = 'Trip ' + CONVERT(VARCHAR(10), @TripLogId) + ' intersecting points with ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName + ': ' + CONVERT(varchar(10), @rc)
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogDebug, @Section

SET @Message = 'Save intersects with ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogInfo, @Section

INSERT INTO @Intersects (LogId, NodeId, ogr_fid)
SELECT 
    LogId,
    NodeId,
    @GeographyId
FROM #Logs
WHERE IsIntersect = 1

What I'm looking for here a way to visualise the zones vs the points so I can tell if they are broken or not. I don't know how to compare fields in the built in spatial results tab.

Answer 1

Below is the final script I used, in case anyone needs the help.

@Shapes is a table variable. It is filled with an insert, and then looped through by a cursor.

Each cursor loop the @Shape is tested for validity. If invalid it has MakeValid() run. It is then tested to see if it is a MULTIPOLYGON . If it is, it is broken into @SubShapes. If it is a single shape, that shape is assigned to the 1st @SubShape.

Each @SubShape is tested for Validation and Area. If required, MakeValid() and ReorientObject() are run.

After all @SubShapes of a @Shape are run, if @Shape is a MULTIPOLYGON it is put back together. This is done via string concatenation and manipulation of the WellKnownText. This is the bit I was least happy with. It feels like a gap in the geography API.

Once all the @Shapes are processed, they are written back to the real table

usp_Log_Error is a an in-house logging procedure. Feel free to comment it out

SET NOCOUNT ON

DECLARE @ProcessName varchar(50) = 'ST_Geography_Validation'
DECLARE @SessionName varchar(50) = 'ST_GV_' + CONVERT(VARCHAR(19), GETDATE(), 126)
DECLARE @Message varchar(1000)

DECLARE @Id int 
DECLARE @Name varchar(254) 
DECLARE @Shape geography
DECLARE @ValidationDetail varchar(200)
DECLARE @ShapeCount int

DECLARE @SubValidationDetail varchar(200)
DECLARE @SubShape geography
Declare @Area float
DECLARE @AreaStr varchar(20)

DECLARE @i int

DECLARE @IsMulti int = 0
DECLARE @MultiWKT varchar(MAX)
DECLARE @SingleWKT varchar(MAX)

DECLARE @Shapes TABLE
(
    Id int,
    Name varchar(254),
    Shape geography
)

INSERT INTO @Shapes (Id, Name, Shape)
SELECT
    ogr_fid,
    name,
    geog4326
FROM safertogether2


DECLARE @SubShapes TABLE
(
    SubId int PRIMARY KEY,
    SubShape geography,
    Processed bit default 0
)

DECLARE _c CURSOR FOR
    SELECT 
        s.Id, 
        s.name, 
        s.Shape,
        0,
        '',
        ''
    FROM @Shapes s

    OPEN _c  
        FETCH NEXT FROM _c INTO @Id, @Name, @Shape, @IsMulti, @MultiWKT, @SingleWKT

        WHILE @@FETCH_STATUS = 0  
        BEGIN

            DELETE FROM @SubShapes

            SET @ValidationDetail = @Shape.IsValidDetailed() -- This is for the base shape
            IF @ValidationDetail <> '24400: Valid'
            BEGIN
                SET @Shape = @Shape.MakeValid()
            END

            SET @ShapeCount = @Shape.STNumGeometries()
            IF @ShapeCount > 1
            BEGIN
                RAISERROR('%i:%s: is MULTIPOLYGON. Break into subshapes.', 0, 1, @Id, @Name) WITH NOWAIT
                SET @IsMulti = 2
                INSERT INTO @SubShapes(SubId, SubShape)
                SELECT 
                    smg.Id,
                    smg.Geog 
                FROM [dbo].uf_SplitMultiGeography(@Shape) smg
            END
            ELSE IF @ShapeCount = 1
            BEGIN 
                RAISERROR('%i:%s: is SINGLE. Insert Shape whole as Id 1.', 0, 1, @Id, @Name) WITH NOWAIT
                SET @IsMulti = 1
                INSERT INTO @SubShapes(SubId, SubShape) 
                VALUES (1, @Shape)
            END
            ELSE
            BEGIN 
                RAISERROR('%i:%s: is NONE.', 0, 1, @Id, @Name) WITH NOWAIT
                SET @IsMulti = 0
            END

            SET @i = 1
            WHILE @i <= @ShapeCount
            BEGIN
                SELECT @SubShape = s.SubShape
                FROM @SubShapes s
                WHERE SubId = @i

                SET @SubValidationDetail = @SubShape.IsValidDetailed()
                RAISERROR('%i:%s:%i: .IsValidDetailed(): %s', 0, 1, @Id, @Name, @i, @SubValidationDetail) WITH NOWAIT

                IF @SubValidationDetail <> '24400: Valid'
                BEGIN
                    RAISERROR('%i:%s:%i: .MakeValid()', 0, 1, @Id, @Name, @i) WITH NOWAIT
                    SET @SubShape = @SubShape.MakeValid() -- This is for the subshape
                END

                BEGIN TRY
                    SET @Area = @SubShape.STArea()
                END TRY
                BEGIN CATCH
                    SET @Message = 'While getting the Area of a subshape'
                    EXEC usp_Log_Error @ProcessName, @SessionName, @Message
                END CATCH
            
                SET @AreaStr = LTRIM(STR(@Area, 20, 5))
                RAISERROR('%i:%s:%i: .STArea(): Pre: %s', 0, 1, @Id, @Name, @i, @AreaStr) WITH NOWAIT
                IF @Area > 510000000000000
                BEGIN
                    RAISERROR('%i:%s:%i: .ReorientObject()', 0, 1, @Id, @Name, @i) WITH NOWAIT
                    SET @SubShape = @SubShape.ReorientObject()
                    SET @Area = @SubShape.STArea()
                    SET @AreaStr = LTRIM(STR(@Area, 20, 5))
                    RAISERROR('%i:%s:%i: .STArea(): Post: %s', 0, 1, @Id, @Name, @i, @AreaStr) WITH NOWAIT
                END

                UPDATE @SubShapes 
                SET SubShape = @SubShape
                WHERE SubId = @Id

                IF @IsMulti = 2
                BEGIN
                    RAISERROR('%i:%s:%i: MULTIPOLYGON member.  Process WKT for recombination', 0, 1, @Id, @Name, @i) WITH NOWAIT
                    BEGIN TRY
                        SET @SingleWKT = @SubShape.STAsText()
                    END TRY
                    BEGIN CATCH
                        SET @Message = 'While getting the WKT of a subshape'
                        EXEC usp_Log_Error @ProcessName, @SessionName, @Message
                    END CATCH

                    SET @SingleWKT = REPLACE(@SingleWKT, 'POLYGON', '')
                    IF (@i <> @ShapeCount)
                    BEGIN
                        SET @SingleWKT = REPLACE(@SingleWKT, '))', ')),')
                    END
                    SET @MultiWKT = @MultiWKT + @SingleWKT
                END
                SET @i = @i + 1
            END

        
            IF @IsMulti = 2
            BEGIN
                RAISERROR('%i:%s: Is a MULTIPOLYGON.  Recreate from updated WKT.', 0, 1, @Id, @Name, @i) WITH NOWAIT
                SET @MultiWKT = 'MULTIPOLYGON(' + @MultiWKT + ')'
                BEGIN TRY
                    SET @Shape = geography::STGeomFromText(@MultiWKT, 4326)
                END TRY
                BEGIN CATCH
                    SET @Message = 'While creating a multipolygon'
                    EXEC usp_Log_Error @ProcessName, @SessionName, @Message
                END CATCH
            END
            ELSE IF @IsMulti = 1
            BEGIN
                RAISERROR('%i:%s: Is SINGLE.  Set Shape to subshape', 0, 1, @Id, @Name, @i) WITH NOWAIT
                SET @Shape = @SubShape
            END
            ELSE
            BEGIN
                RAISERROR('%i:%s: Is NONE.  Set Shape to none', 0, 1, @Id, @Name, @i) WITH NOWAIT
            END

            RAISERROR('%i:%s: Update temp table', 0, 1, @Id, @Name) WITH NOWAIT
            UPDATE @Shapes
            SET Shape = @Shape
            WHERE Id = @Id

            FETCH NEXT FROM _c INTO @Id, @Name, @Shape, @IsMulti, @MultiWKT, @SingleWKT
        END
    CLOSE _c
DEALLOCATE _c

select 
    st.ogr_fid, 
    st.name, 
    st.geog4326 as ShapePre, 
    s.Shape as ShapePost,
    CASE 
        WHEN  st.geog4326.IsValidDetailed() = '24400: Valid' then st.geog4326.STArea()
        ELSE st.geog4326.MakeValid().STArea()
    END as AreaPre,
    CASE 
        WHEN  s.Shape.IsValidDetailed() = '24400: Valid' then s.Shape.STArea()
        ELSE s.Shape.MakeValid().STArea()
    END as AreaPost  
from safertogether2 st
join @Shapes s ON s.Id = st.ogr_fid

RAISERROR('Update back to safertogether table', 0, 1) WITH NOWAIT
UPDATE st
SET geog4326 = s.Shape
FROM @Shapes s
JOIN safertogether2 st ON s.Id = st.ogr_fid


RAISERROR('End of Script', 0, 1, @Id, @Name) WITH NOWAIT

Answer 2

It is likely the problem is in conversion of shapefile (planar map) data to SQL Server Geography (spherical) model. These are easy to do wrong and get inverted (complimentary) polygons, i.e. polygons describing all Earth area except the intended shape.

Also see couple similar cases, they were with Google BigQuery, which uses similar semantics of spherical oriented polygons:

• Query times out after 6 hours, how to optimize it?
• How to get the correct centroid of a bigquery polygon with st_centroid

You can check if the polygons were ingested correctly to the database by computing their area - it will be huge (order of 5e14 m^2) for inverted polygons, whereas you should get reasonable result for correctly ingested ones.

SQL Server has ReorientObject() method that reverses the polygon orientation. Be careful with it though, as orientation interacts with validity - so it is better to use ReorientObject() before MakeValid(). Best of all of course is to fix the ingestion method.

Alternatively, you can use SQL Server Geometry (rather than Geography) type that uses planar map semantics. Working with planar data has some issues, but assuming you don't cross anti-meridian, and have basic understanding of which projections are good for which kind of computation, you should typically be OK.

Or switch to a database that uses spherical mode, but can interpret and convert planar data correctly. E.g. BigQuery understands that GeoJson geometries are planar, and ignores their orientation when converting them to the internal spherical Geography model, so you never have this problem if you load Geometry as GeoJson fields.