SQL Geography复杂形状内切矩形计算与地理边界重叠分析

Hey there, let's work through these two SQL Geography tasks step by step—both are super relevant when dealing with spatial data in SQL Server, so I'll share practical, actionable solutions for each.

SQL Server doesn’t have a native function to compute this directly, so we’ll use an iterative grid-based approach to approximate the largest rectangle that fits entirely inside your target shape. Here’s how to pull it off:

Approach Overview #

1. Start with the bounding box of your geography shape (using STEnvelope()) to define our initial search area.
2. Iteratively divide this area into smaller grid cells, checking which cells are fully contained within the original shape.
3. Track the largest valid cell (by area) and refine the search around it to get a more precise result.

Example Code #

DECLARE @TargetShape GEOGRAPHY = (SELECT Geog FROM YourShapeTable WHERE ShapeID = 1);
DECLARE @MinPrecision FLOAT = 100; -- Minimum cell size in meters (adjust based on your needs)
DECLARE @CurrentGridSize FLOAT = @TargetShape.STEnvelope().STLength() / 10; -- Start with 10x10 grid

-- Variables to track the best rectangle found
DECLARE @BestRectangle GEOGRAPHY;
DECLARE @BestArea FLOAT = 0;

WHILE @CurrentGridSize >= @MinPrecision
BEGIN
    -- Generate grid cells within the bounding box
    WITH Grid AS (
        SELECT 
            x, y,
            @TargetShape.STEnvelope().STPointN(1).STX + (x * @CurrentGridSize) AS X1,
            @TargetShape.STEnvelope().STPointN(1).STY + (y * @CurrentGridSize) AS Y1,
            @TargetShape.STEnvelope().STPointN(1).STX + ((x+1) * @CurrentGridSize) AS X2,
            @TargetShape.STEnvelope().STPointN(1).STY + ((y+1) * @CurrentGridSize) AS Y2
        FROM (SELECT TOP 1000 ROW_NUMBER() OVER (ORDER BY (SELECT NULL)) -1 AS x FROM sys.all_columns) x
        CROSS JOIN (SELECT TOP 1000 ROW_NUMBER() OVER (ORDER BY (SELECT NULL)) -1 AS y FROM sys.all_columns) y
        WHERE 
            @TargetShape.STEnvelope().STPointN(1).STX + (x * @CurrentGridSize) <= @TargetShape.STEnvelope().STPointN(3).STX
            AND @TargetShape.STEnvelope().STPointN(1).STY + (y * @CurrentGridSize) <= @TargetShape.STEnvelope().STPointN(3).STY
    ),
    CandidateRectangles AS (
        SELECT 
            GEOGRAPHY::STPolyFromText(
                'POLYGON((' + 
                CAST(X1 AS VARCHAR) + ' ' + CAST(Y1 AS VARCHAR) + ', ' +
                CAST(X2 AS VARCHAR) + ' ' + CAST(Y1 AS VARCHAR) + ', ' +
                CAST(X2 AS VARCHAR) + ' ' + CAST(Y2 AS VARCHAR) + ', ' +
                CAST(X1 AS VARCHAR) + ' ' + CAST(Y2 AS VARCHAR) + ', ' +
                CAST(X1 AS VARCHAR) + ' ' + CAST(Y1 AS VARCHAR) + '))',
                @TargetShape.STSrid
            ) AS Rect,
            (X2 - X1) * (Y2 - Y1) AS Area
        FROM Grid
    )
    SELECT TOP 1 @BestRectangle = Rect, @BestArea = Area
    FROM CandidateRectangles
    WHERE Rect.STWithin(@TargetShape) = 1
    ORDER BY Area DESC;

    -- Reduce grid size for next iteration
    SET @CurrentGridSize = @CurrentGridSize / 2;
END

-- Output the result
SELECT @BestRectangle AS MaxInscribedRectangle, @BestArea AS RectangleArea;

Notes #

• Adjust @MinPrecision based on how accurate you need the rectangle to be (smaller values = more precise but slower).
• For extremely complex shapes, you might want to add a step to first compute the convex hull (STConvexHull()) to narrow down the search area.

This is a straightforward spatial join task once you have your DMA and zip code geography data aligned (make sure both use the same SRID, usually 4326 for WGS84). We’ll break this into three parts: fully contained zips, partially overlapping zips, and overlap percentage.

Prerequisites #

• Two tables: DMA_Table (with DMA_ID, DMA_Name, Geog columns) and ZipCode_Table (with ZipCode, Geog columns).
• Spatial indexes on both Geog columns to speed up queries (critical for large datasets).

1. Find Zip Codes Fully Contained in a DMA #

SELECT
    d.DMA_Name,
    z.ZipCode
FROM DMA_Table d
INNER JOIN ZipCode_Table z
    ON z.Geog.STWithin(d.Geog) = 1
ORDER BY d.DMA_Name, z.ZipCode;

2. Find Partially Overlapping Zip Codes + Overlap Percentage #

SELECT
    d.DMA_Name,
    z.ZipCode,
    -- Calculate overlap area (in square meters by default; convert if needed)
    z.Geog.STIntersection(d.Geog).STArea() AS OverlapArea_SqMeters,
    -- Calculate percentage of the zip code that lies within the DMA
    ROUND(
        (z.Geog.STIntersection(d.Geog).STArea() / z.Geog.STArea()) * 100,
        2
    ) AS OverlapPercentage
FROM DMA_Table d
INNER JOIN ZipCode_Table z
    ON z.Geog.STIntersects(d.Geog) = 1
WHERE z.Geog.STWithin(d.Geog) = 0 -- Exclude fully contained zips
ORDER BY d.DMA_Name, OverlapPercentage DESC;

Performance Tips #

• Add spatial indexes: CREATE SPATIAL INDEX SIX_DMA_Geog ON DMA_Table(Geog); and same for the zip code table.
• Use STIntersects first to filter down the dataset before running more expensive STWithin or STIntersection calculations.
• If you need area in square miles, divide the result by 2589988.110336 (since 1 square mile = ~2.59 million square meters).

内容的提问来源于stack exchange，提问作者Jake McGraw