SQL Server STWithin, STIntersects inaccurate

Here's my setup. I define a polygon and a point close to the middle of that polygon. Using STWithin or STIntersects nothing is returned.

create table #gtest ([Lat] [float], [Lng] [float], [GeoPoint]  AS ([geography]::Point([Lat],[Lng],(4326))))
insert #gtest(Lat, Lng) values (36.3893, -102.796)

declare @swLng float, @swLat floa, @neLng float, @neLat float
declare @g geography, @poly varchar(max)

select @swLng=-147.45514585847684, @swLat=28.841195680455485, @neLng=-59.86599816500535, @neLat=56.14838208179307
set @poly = 'POLYGON((' + cast(@swLng as varchar) + ' ' +  cast(@swLat as varchar) + ',' + cast(@neLng as varchar) + ' ' + cast(@swLat as varchar) + ',' + cast(@neLng as varchar) + ' ' + cast(@neLat as varchar) + ',' + cast(@swLng as varchar) + ' ' + cast(@neLat as varchar) + ',' +  cast(@swLng as varchar) + ' ' +  cast(@swLat as varchar) + '))'
set @g = geography::STGeomFromText(@poly, 4326);

select * from #gtest where GeoPoint.STWithin(@g) = 1 -- no results
select * from #gtest where GeoPoint.STIntersects(@g) = 1 -- no results either

But if I decrease the polygon size, the point starts to fit in it from sql server perspective.

select @swLng=-134.35054327724572, @swLat=31.64106010841097, @neLng=-71.77427705012154, @neLat=51.66708602548192

UPD: From the comments I can see that the polygon is a section of a sphere, which makes sense but how do I deal with Mercator projection? I'm plotting some spots using Mapbox and that's what I'm getting, using Mapbox's boundaries:

Solution

If you want the boundaries to follow lines of latitude or otherwise follow straight lines on a Mercator projection, I believe you will need to approximate it (to an acceptable tolerance) with a series of line segments by adding intermediate points between each pair or existing vertices.

For this discussion, I will use (x, y) coordinates notation instead of (Longitude, Latitude). The approach below is also an approximation that assumes that the error to be corrected is equivalent to the difference between an arc and a chord in Euclidean space.

If your edge start and end coordinates are defined as {(xs,ys),(xe,ye)}, you can define a Euclidean midpoint as (xm,ym), where xm = (xs+xe)/2 and xm = (xs+xe)/2. You can then use STDistance to calculate the distance h between that point and the line (great-circle arc) defined by your original endpoints. (See correction below.)

For a given arc of radius R and angle θ, the chord height h is equal to R(1-cos(θ/2)), where θ is the angle of the arc. For small arcs, this height is roughly proportional to the angle-squared, so to achieve a tolerance t (in meters), you can calculate N ≈ SQRT(h/t). Because this is an approximation, you might need to add a fudge-factor and calculate something like N = CEIL(1.1 * h / t).

From there you can calculate the intermediate points (xi,yi) for i in 1..N-1, where xi = ((N-i) * xs + i * xe) / N and yi = ((N-i) * ys + i * ye) / N. You would insert these points between y our original start and end points. (See correction below.)

For boundary segments that extend north/south (constant longitude) and others where the original h is less than the tolerance t, you may find that the calculated N is 0 or 1, in which case you do not need to insert any intermediate points.

Correction: The ym (midpoint latitude) and yi (interpolated latitude) calculations above did not account for the vertical stretching of the projected Mercator coordinates. To properly handle this, the start/end latitude values can be mapped to the projected vertical position, calculation of the intermediate points can be performed in in that projected space, and then a reverse projection applied to the resulting coordinates. The x/longitude components are not affected. The calculations would then look something like:

ysp = project(ys)
yep = project(ye)
ymp = (ysp + yep) / 2
ym = unproject(ymp)
yip = ((N-i) * ysp + i * yep) / N
yi = unproject(yip) For horizontal or vertical boundaries (constant y or x), this updated calculation produces similar results, but for an angled edge, the original calculation might produce a curved diagonal path.

See Mercator projection transformations

Also, longitude calculations need to recognize and handle +/- 180 degree longitude crossover situations - first by recognizing when the start/end distance is shorter across the 180 degree line, and performing intermediate point calculations using some variation of modulo arithmetic.

Everything I have written above is untested at this point, but I am working on some sample code. I really went down the rabbit hole, It was an interesting exercise. Here is what I have so far.

First a few helper functions for mapping latitude values to and from the projected Mercator vertical axis and for interpolating latitude and longitude values. Latitude interpolation is performed in the projected vertical axis. Longitude interpolation include 180 degree crossover handling.

CREATE FUNCTION LatToMercatorProjectedLat(@Lat FLOAT)
RETURNS FLOAT
AS BEGIN
    -- Degenerate case @Lat ~ +/- 90 degrees
    IF ABS(@Lat) > 89.9999999999
        RETURN 9999 * SIGN(@Lat)

    -- https://en.wikipedia.org/wiki/Mercator_projection#Derivation
    -- Modified to include degrees/radians conversion
    RETURN LOG(TAN((45 + 0.5 * @Lat) * PI() / 180)) / PI() * 180
END

CREATE FUNCTION MercatorProjectedToDegreesLat(@Projected FLOAT)
RETURNS FLOAT
AS BEGIN
    -- Degenerate case @Lat ~ +/- 90 degrees (round trip)
    IF ABS(@Projected) >= 9999
        RETURN 90 * SIGN(@Projected)

    -- https://en.wikipedia.org/wiki/Mercator_projection#Derivation
    -- Modified to include degrees/radians conversion
    RETURN (ATAN(EXP(@Projected * PI() / 180)) / PI() * 180 - 45) * 2
END

CREATE FUNCTION InterpolateLatitude(
    @Latitude1 FLOAT,  -- Start
    @Latitude2 FLOAT,  -- End
    @I         FLOAT,  -- Interpolation factor 0 (start) to N (end)
    @N         FLOAT,  -- Number of steps
    @Mercator  BIT     -- If 1, adjust for Mercator projection
)
RETURNS FLOAT -- Interpolated resultS
BEGIN
    IF @Mercator = 1
    BEGIN
        SET @Latitude1 = dbo.LatToMercatorProjectedLat(@Latitude1)
        SET @Latitude2 = dbo.LatToMercatorProjectedLat(@Latitude2)  
    END

    DECLARE @Fraction FLOAT = @I / NULLIF(@N, 0)
    DECLARE @Result FLOAT = @Latitude1 * (1 - @Fraction) + @Latitude2 * @Fraction

    IF @Mercator = 1
    BEGIN
        SET @Result = dbo.MercatorProjectedToDegreesLat(@Result)
    END

    RETURN @Result
END

CREATE FUNCTION InterpolateLongitude(
    @Longitude1 FLOAT,  -- Start longitude in degrees
    @Longitude2 FLOAT,  -- End longitude in degrees
    @I          FLOAT,  -- Interpolation factor 0 (start) to N (end)
    @N          FLOAT   -- Number of steps
)
RETURNS FLOAT -- Interpolated result, adjusted to be in +/- 180 range.
AS
BEGIN
    -- For longitudes that cross the +/-180 mark, adjust
    -- @LatOrLong2 to be within +/- 180 degrees from @LatOrLong1.
    DECLARE @AdjustedLongitude2 FLOAT = CASE
        WHEN @Longitude2 < @Longitude1 - 180 THEN @Longitude2 + 360
        WHEN @Longitude2 > @Longitude1 + 180 THEN @Longitude2 - 360
        ELSE @Longitude2
        END
    -- Calculate interpolated latitude or longitude and adjust to be in +/- 180 range.
    DECLARE @Fraction FLOAT = @I / NULLIF(@N, 0)
    DECLARE @Temp FLOAT =
        @Longitude1 * (1 - @Fraction)
        + @AdjustedLongitude2 * @Fraction
        + 180
    DECLARE @NormalizedResult FLOAT =
        @Temp
        - 360 * FLOOR(@Temp / 360)
        - 180
    RETURN @NormalizedResult
END

Now for the main event. The following will take an existing polygon geometry object and:

Break it down into component boundaries, vertices, and edges.
Calculate the midpoint, estimated error, and number of required midpoints for each edge.
Generate the required number of additional edge midpoints.
Combine all the data back into a new polygon geometry object.

CREATE FUNCTION AdjustPolygonEdgesForMercadorProjection(
    @Polygon Geography,
    @Tolerance FLOAT)
RETURNS Geography
AS
BEGIN
    DECLARE @Result Geography
    DECLARE @Fudge FLOAT   = 1.1  -- Allowance for approximations
    DECLARE @Mercator BIT  = 1    -- If 1, calc using Mercator projected latitude

    DECLARE @SRID INT = @Polygon.STSrid -- Typically 4326 for WGS84

    ;WITH Curves AS (
        -- A polygon will have one outer boundary curve plus zero or more
        -- inner boundary curves (holes)
        SELECT 
            S.value AS CurveNum,
            @Polygon.RingN(S.value) AS Curve
        FROM GENERATE_SERIES(1, @Polygon.NumRings()) S
    ),
    Vertices AS (
        -- Each curve (outer and inner boundary) should have 4 or more vertices
        -- (last = first) representing 3 or more edges.
        SELECT
            C.CurveNum,
            S.value AS PointNum,
            C.Curve.STPointN(S.value).Long AS Lng,
            C.Curve.STPointN(S.value).Lat AS Lat
        FROM Curves C
        CROSS APPLY GENERATE_SERIES(1, C.Curve.STNumPoints()) S
    ),
    Edges AS (
        -- Extract edges.
        SELECT
            V1.CurveNum,
            V1.PointNum AS EdgeNum,
            V1.Lng AS Lng1, V1.Lat AS Lat1,
            V2.Lng AS Lng2, V2.Lat AS Lat2
        FROM Vertices V1
        JOIN Vertices V2
            ON V2.CurveNum = V1.CurveNum
            AND V2.PointNum = V1.PointNum + 1
    ),
    InterpolatedPoints AS (
        -- Original points, including last which equals first
        SELECT V.CurveNum, V.PointNum, 0 AS InterNum, V.Lng, V.Lat
        FROM Vertices V

        UNION ALL

        -- Interpolated points (if needed)
        SELECT E.CurveNum, E.EdgeNum AS PointNum, I.InterNum, I.LngI, I.LatI
        FROM Edges E
        CROSS APPLY (
            -- Calculate midpoint (*not* great-circle)
            SELECT
                dbo.InterpolateLongitude(E.Lng1, E.Lng2, 1, 2) AS LngM,
                dbo.InterpolateLatitude(E.Lat1, E.Lat2, 1, 2, @Mercator) AS LatM
        ) M
        CROSS APPLY (
            -- Define the great circle line (spheroid arc) text
            SELECT
                CONCAT(
                    'LINESTRING(',
                    CONVERT(VARCHAR, E.Lng1, 3), ' ', CONVERT(VARCHAR, E.Lat1, 3),
                    ', ',
                    CONVERT(VARCHAR, E.Lng2, 3), ' ', CONVERT(VARCHAR, E.Lat2, 3),
                    ')'
                    ) AS LineText
        ) GT
        CROSS APPLY (
            -- Define the great circle line (arc) and mercator midpoint geometry objects
            SELECT
                Geography::STLineFromText(GT.LineText, @SRID) AS Line,
                Geography::Point(M.LatM, M.LngM, @SRID) AS MidPoint
        ) G
        CROSS APPLY (
            -- Calculate the distance from the mercator midpoint to the
            -- great-circle line
            SELECT G.Line.STDistance(G.MidPoint) AS Distance
        ) D
        CROSS APPLY (
            -- Calculate the number of segments needed to reduce error to tolerance
            SELECT CAST(CEILING(SQRT(@Fudge * D.Distance / @Tolerance)) AS INT) AS N
        ) N
        CROSS APPLY (
            -- Generate intermediate points (excludes original start/end points)
            SELECT
                S.value AS InterNum,
                dbo.InterpolateLongitude(E.Lng1, E.Lng2, S.value, N.N) AS LngI,
                dbo.InterpolateLatitude(E.Lat1, E.Lat2, S.value, N.N, @Mercator) AS LatI
            FROM GENERATE_SERIES(1, N.N - 1) S
        ) I
        WHERE N.N >= 2
    ),
    -- We have what we need, it is time to reassemble the results
    PointText AS (
        -- Format text for each point in the form "X Y", preserving max precision
        SELECT
            IP.CurveNum, IP.PointNum, IP.InterNum,
            CAST(CONCAT(
                CONVERT(VARCHAR, IP.Lng, 3),
                ' ',
                CONVERT(VARCHAR, IP.Lat, 3)
            ) AS VARCHAR(MAX)) AS PointText
        FROM InterpolatedPoints IP
    ),
    CurveText AS (
        -- Combine points into curves of the form "(X1 Y1, ..., Xn Yn)"
        -- one row for outer polygon boundary (CurveNum = 1)
        -- plus zero or more rows for inner boundaries (CurveNum >= 2)
        SELECT
            PT.CurveNum,
            CONCAT(
                '(',
                STRING_AGG(PT.PointText, ', ')
                    WITHIN GROUP(ORDER BY PT.PointNum, PT.InterNum),
                ')'
            ) AS CurveText
        FROM PointText PT
        GROUP BY PT.CurveNum
    ),
    GeoText AS (
        -- Combine outer and (optional) inner boundaries into polygon text
        -- "POLYGON((X1 Y1, ..., Xn Yn))" or "POLYGON((...), (...), ...)"
        SELECT
            CONCAT(
                'POLYGON(',
                STRING_AGG(CT.CurveText, ', ')
                     WITHIN GROUP(ORDER BY CT.CurveNum),
                ')'
            ) AS PolygonText
        FROM CurveText CT
    ),
    Geo AS (
        -- Convert polygon text to polygon
        SELECT Geography::STGeomFromText(GT.PolygonText, @SRID) AS Polygon
        FROM GeoText GT
    )
    SELECT @Result = G.Polygon
    FROM Geo G

    RETURN @Result
END

I have used CTEs and CROSS APPLYs extensively to incrementally build up intermediate results and to avoid excessive nested subselects. It may be wordy, but my intent is to try to make it readable by implementing the logic in small steps.

At this point, it appears to work, but I have only performed limited testing. If used, an more thorough test should be performed.

See this db<>fiddle for a demo and a few extra tests.