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

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

@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

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

@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.

-- 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

@Polygon Geography,

@Tolerance FLOAT)

RETURNS Geography

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,

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 float, @neLng float, @neLat float

select @swLng=-147.45514585847684, @swLat=28.841195680455485, @neLng=-59.86599816500535, @neLat=56.14838208179307

declare @g geography, @poly varchar(max)

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 @g

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

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

declare @tolerance float = 100.0 -- 100 meters tolerance

set @g = dbo.AdjustPolygonEdgesForMercadorProjection(@g, @tolerance)

select @g

select * from #gtest where GeoPoint.STWithin(@g) = 1 -- Match!

select * from #gtest where GeoPoint.STIntersects(@g) = 1 -- Also a match!

L.Lat,

dbo.LatToMercatorProjectedLat(L.Lat) AS MercatorProjected,

dbo.MercatorProjectedToDegreesLat(

dbo.LatToMercatorProjectedLat(L.Lat)

) AS Inverse

FROM (

SELECT S.value AS Lat

FROM GENERATE_SERIES(-100, 100, 10) S -- Includes some out-of range values

) L

ORDER BY L.Lat DESC

Lat,

dbo.LatToMercatorProjectedLat(Lat) AS MercatorProjected,

dbo.MercatorProjectedToDegreesLat(

dbo.LatToMercatorProjectedLat(Lat)

) AS Inverse

FROM (

VALUES

(89), (89.99), (89.9999), (89.999999), (89.99999999), (89.9999999999),

(90)

) V(Val)

CROSS APPLY (VALUES (V.Val), (-V.Val)) L(Lat)

ORDER BY Lat DESC

*,

MercatorProjectedLat - LAG(MercatorProjectedLat)

OVER(PARTITION BY TestCaseId ORDER BY I)

AS MercatorProjectedDiff

FROM (

VALUES

(1, 60, 75, 0),

(2, -60, -75, 1),

(3, 60, 75, 2),

(4, -60, -75, 3),

(5, -10, 80, 4), -- 0 degree crossover

(6, -120, 150, 5) -- 180 longitude crossover (latititude is part out of range)

) V(TestCaseId, L1, L2, N)

OUTER APPLY (

SELECT

S.value AS I,

dbo.InterpolateLongitude(L1, L2, S.value, N) Long,

dbo.InterpolateLatitude(L1, L2, S.value, N, 0) EuclideanLat,

dbo.InterpolateLatitude(L1, L2, S.value, N, 1) MercatorLat

FROM GENERATE_SERIES(0, N) S

) Calc

OUTER APPLY (

SELECT

dbo.LatToMercatorProjectedLat(MercatorLat) MercatorProjectedLat

) Calc2

DECLARE @Tolerance FLOAT = 1000 -- 1000 meters (1km)

DECLARE @SRID INT = 4326 -- WGS84

DECLARE @N INT = 1000 -- Points to test per edge

G.Polygon.STNumPoints() AS NPoints,

G2.Polygon2.STNumPoints() AS NPoints2,

G3.Polygon3.STNumPoints() AS NPoints3,

--G.Polygon.STArea() * 1e-6 AS SqKm1,

--G2.Polygon2.STArea() * 1e-6 AS SqKm2,

MD.*,

G2.Polygon2

FROM (

VALUES

(-147.45514585847684, 28.841195680455485, -59.86599816500535, 56.14838208179307),

(-134.35054327724572, 31.64106010841097, -71.77427705012154, 51.66708602548192),

(170, 28, -59, 56)

) V(Lng1, Lat1, Lng2, Lat2)

CROSS APPLY (

SELECT

CONVERT(VARCHAR, Lng1, 3) AS Lng1Text,

CONVERT(VARCHAR, Lat1, 3) AS Lat1Text,

CONVERT(VARCHAR, Lng2, 3) AS Lng2Text,

CONVERT(VARCHAR, Lat2, 3) AS Lat2Text

) T

CROSS APPLY (

SELECT CONCAT(

'POLYGON((',

Lng1Text, ' ', Lat1Text, ',',

Lng2Text, ' ', Lat1Text, ',',

Lng2Text, ' ', Lat2Text, ',',

Lng1Text, ' ', Lat2Text, ',',

Lng1Text, ' ', Lat1Text,

'))'