STCentroid*: Alternate Functions for Compute a Centroid
The top download from my website is for my package of centroid algorithms for Oracle sdo_geometry.
There is nothing wrong with the STCentroid function implemented in SQL Server Spatial (for one thing it handles geography objects natively), but it seems that having a choice of centroid algorithms is a big thing for users of spatial data within databases.
The original algorithm for getting the centroid of a polygon involves use of LAG/LEAD SQL99 analytic functions that have not been a part of SQL Server until Denali.
For the polygon centroid function (STCentroid_P) one can ask the function to generate a centroid:
- In the middle of the MBR of the geometry, or
- In a weighted position based on coordinate density, or
- With a defined or fixed (supplied) “seed” X ordinate position.
Here are the functions.
Centroid for MultiPoint Geometries
SQL Server Spatial returns null when executing STCentroid against multipoint geometries, so let’s build our own.
USE [GISDB] -- Change to your database GO -- MultiPoint Centroid Creator -- DROP FUNCTION [dbo].[STCENTROID_P] ; CREATE FUNCTION [dbo].[STCENTROID_P] ( @p_geometry GEOMETRY, @p_round_x INT = 3, @p_round_y INT = 3, @p_round_z INT = 3 ) RETURNS GEOMETRY AS BEGIN DECLARE @v_x FLOAT = 0.0, @v_y FLOAT = 0.0, @v_z FLOAT, @v_geomn INT, @v_part_geom geometry, @v_geometry geometry, @v_gtype nvarchar(MAX); BEGIN IF ( @p_geometry IS NULL ) RETURN @p_geometry; SET @v_gtype = @p_geometry.STGeometryType(); IF ( @v_gtype = 'MultiPoint' ) BEGIN -- Get parts of multi-point geometry -- SET @v_geomn = 1; WHILE ( @v_geomn <= @p_geometry.STNumGeometries() ) BEGIN SET @v_part_geom = @p_geometry.STGeometryN(@v_geomn); SET @v_x = @v_x + @v_part_geom.STX SET @v_y = @v_y + @v_part_geom.STY SET @v_geomn = @v_geomn + 1; END; SET @v_geometry = geometry::Point(@v_x/@p_geometry.STNumGeometries(), @v_y/@p_geometry.STNumGeometries(), @p_geometry.STSrid); END ELSE SET @v_geometry = @p_geometry; RETURN @v_geometry; END; END GO
Now let’s test noting that SQL Server’s STCentroid always returns NULL for MultiPoint geometries
WITH multiPoint AS (
SELECT geometry::STGeomFromText('MULTIPOINT((0 0),(100 0),(100 100),(0 100),(150 110),(150 150),(110 150),(110 110))',0) AS geom
)
SELECT 'O' AS label, geom.STBuffer(2) AS geom FROM multiPoint
UNION ALL
SELECT 'M' AS label, dbo.STCentroid_P(geom,3,3,3).STBuffer(5) AS geom FROM multiPoint
UNION ALL
SELECT 'S' AS label, geom.STCentroid().STBuffer(5) AS geom FROM multiPoint;
Which looks like this:
Centroid for Area Geometries
DROP FUNCTION [dbo].[STCENTROID_A];
CREATE FUNCTION [dbo].[STCENTROID_A] (
@p_geometry geometry,
@p_start INT = 0, /* 0 = use average of all Area's vertices for starting X centroid calculation
1 = centre X of MBR
2 = User supplied starting seed X */
@p_seed_x FLOAT = NULL,
@p_round_x INT = 3,
@p_round_y INT = 3,
@p_round_z INT = 3
)
RETURNS geometry
AS
BEGIN
DECLARE
@v_cx FLOAT,
@v_cy FLOAT,
@v_seed geometry,
@v_wkt nvarchar(MAX);
BEGIN
IF ( @p_geometry IS NULL )
RETURN @p_geometry;
IF ( @p_start NOT IN (0,1,2) )
RETURN CAST('Starting position must be 0, 1 or 2.' AS VARCHAR(MAX)); -- geometry);
IF ( @p_start = 2 AND @p_seed_x IS NULL )
RETURN CAST('Seed X value not provided.' AS VARCHAR(MAX)); -- geometry);
-- Create starting seed point geometry
--
SET @v_seed = CASE @p_start
WHEN 0 THEN @p_geometry
WHEN 1 THEN geometry::Point(@p_geometry.STEnvelope().STPointN(1).STX,
@p_geometry.STEnvelope().STPointN(1).STY,
@p_geometry.STSrid).STUnion(
geometry::Point(@p_geometry.STEnvelope().STPointN(3).STX,
@p_geometry.STEnvelope().STPointN(3).STY,
@p_geometry.STSrid))
WHEN 2 THEN geometry::Point(@p_seed_x,0,@p_geometry.STSrid)
ELSE @p_geometry
END;
-- Check user seed X between MBR of object
--
IF ( @p_start = 2 AND
( @p_seed_x <= @p_geometry.STEnvelope().STPointN(1).STX OR @p_seed_x >= @p_geometry.STEnvelope().STPointN(3).STX ) )
RETURN CAST('Seed X value not between provided geometry''s MBR.' AS VARCHAR(MAX)); -- geometry);
SELECT TOP 1
@v_cx = cx,
@v_cy = cy
FROM (SELECT z.x AS cx,
z.y + ( ydiff / 2 ) AS cy,
ydiff
FROM (SELECT w.id,
w.x,
w.y,
CASE WHEN w.ydiff IS NULL THEN 0 ELSE w.ydiff END AS ydiff,
CASE WHEN w.id = 1
THEN CASE WHEN w.INOUT = 1
THEN 'INSIDE'
ELSE 'OUTSIDE'
END
WHEN SUM(w.INOUT) OVER (ORDER BY w.id) % 2 = 1
/* Need to look at previous result as inside/outside is a binary switch */
THEN 'INSIDE'
ELSE 'OUTSIDE'
END AS INOUT
FROM (SELECT ROW_NUMBER() OVER (ORDER BY u.y ASC) AS id,
u.x,
u.y,
CASE WHEN u.touchCross IN (0) /* Cross */ THEN 1
WHEN u.touchCross IN (-1,-2,1,2) /* Touch */ THEN 0
ELSE 0
END AS INOUT,
ABS(LEAD(u.y,1) OVER(ORDER BY u.y) - u.y) AS YDiff
FROM (SELECT s.x,
s.y,
/* In cases where polygons have boundaries/holes that touch at a point we need to count them more than once */
CASE WHEN COUNT(*) > 2 THEN 1 ELSE SUM(s.touchcross) END AS touchcross
FROM (SELECT t.x,
t.y,
t.touchcross
FROM (SELECT r.x,
CASE WHEN (r.endx = r.startx)
THEN (r.starty + r.endy ) / 2
ELSE round(r.starty + ( (r.endy-r.starty)/(r.endx-r.startx) ) * (r.x-r.startx),@p_round_y)
END AS y,
CASE WHEN ( r.x = r.startx AND r.x = r.endx )
THEN 99 /* Line is Vertical */
WHEN ( ( r.x = r.startx AND r.x > r.endx )
OR ( r.x = r.endX AND r.x > r.startX ) )
THEN -1 /* Left Touch */
WHEN ( ( r.x = r.endX AND r.x < r.startX )
OR ( r.x = r.startX AND r.x < r.endX ) )
THEN 1 /* Right Touch */
ELSE 0 /* cross */
END AS TouchCross
FROM (SELECT c.x,
round(v.sx,@p_round_x) AS startX,
round(v.sy,@p_round_y) AS startY,
round(v.ex,@p_round_x) AS endX,
round(v.ey,@p_round_y) AS endY
FROM (SELECT round(avg(p.x),@p_round_x) AS x
FROM [dbo].[STDUMPPOINTS](@v_seed) p
) c
CROSS apply
[dbo].[STVECTORIZE](@p_geometry) v
) r
WHERE r.x BETWEEN r.StartX AND r.endx
OR r.x BETWEEN r.endx AND r.startx
) t
) s
GROUP BY s.x,s.y
) u
) w
) z
WHERE z.INOUT = 'INSIDE'
) f
ORDER BY f.ydiff DESC;
RETURN geometry::Point(@v_cx, @v_cy, @p_geometry.STSrid);
END;
END
GO
Now let’s test this function weighting for area (vertex density), MBR (will be the same as the geometry is balanced around its centroid spine of X=2300 ) and via a supplied seed. Comparison to STCentroid is provided.
WITH poly AS (
SELECT geometry::STGeomFromText('POLYGON((2300 -700, 2800 -300, 2300 700, 2800 1100, 2300 1100, 1800 1100, 2300 400, 2300 200, 2100 100, 2500 100, 2300 -200, 1800 -300, 2300 -500, 2200 -400, 2400 -400, 2300 -700),
(2300 1000, 2400 900, 2200 900, 2300 1000))',0) AS geom
)
SELECT 'O' AS Label, geom FROM poly
UNION ALL
SELECT 'A' AS Label, [GISDB].[dbo].[STCentroid_A](geom,0,NULL,3,3,3).STBuffer(15) AS geom FROM poly
UNION ALL
SELECT 'M' AS Label, [GISDB].[dbo].[STCentroid_A](geom,1,NULL,3,3,3).STBuffer(15) AS geom FROM poly
UNION ALL
SELECT 'U' AS Label, [GISDB].[dbo].[STCentroid_A](geom,2,2050,3,3,3).STBuffer(15) AS geom FROM poly
UNION ALL
SELECT 'S' AS Label, geom.STCentroid().STBuffer(15) AS geom FROM poly;
This looks like:
Testing a different geometry shows all four positions more clearly.
WITH weightedPoly AS (
SELECT geometry::STGeomFromText('POLYGON((258.72254365233152 770.97400259630615, 268.79365642517564 739.08214548229967, 278.86476919801976 707.1902883682933, 332.57737065318844 693.76213800450114, 366.14774656266889 676.97695004976094, 426.57442319973364 697.11917559544918, 520.57147574627891 737.40362668682576, 631.35371624756431 744.11770186872184, 829.41893411349884 797.83030332389046, 1547.8249785763801 791.11622814199438, 1205.4071442996797 895.18439346138371, 832.77597170444687 1039.5370098721496, 490.3581374277465 1086.5355361454222, 416.50331042688953 1076.464423372578, 381.25441572193506 1059.6792354178378, 346.00552101698065 1042.8940474630976, 320.82773908487036 1019.3947843264614, 295.64995715276001 995.89552118982499, 287.25736317538986 964.00366407581862, 278.86476919801976 932.11180696181225, 282.2218067889678 891.82735587043567, 277.18625040254574 858.25697996095528, 272.15069401612368 824.68660405147489, 258.72254365233152 770.97400259630615))',0) AS geom
)
SELECT 'O' AS Label, geom FROM weightedPoly
UNION ALL
SELECT 'A' AS Label, [GISDB].[dbo].[STCentroid_A](geom,0,NULL,2,2,2).STBuffer(10) AS geom FROM weightedPoly
UNION ALL
SELECT 'M' AS Label, [GISDB].[dbo].[STCentroid_A](geom,1,NULL,2,2,2).STBuffer(10) AS geom FROM weightedPoly
UNION ALL
SELECT 'U' AS Label, [GISDB].[dbo].[STCentroid_A](geom,2,1200,2,2,2).STBuffer(10) AS geom FROM weightedPoly
UNION ALL
SELECT 'S' AS Label, geom.STCentroid().STBuffer(10) AS geom FROM weightedPoly;
Looks like this:
Centroid for Linear geometries
DROP FUNCTION [dbo].[STCentroid_L];
CREATE FUNCTION [dbo].[STCENTROID_L] (
@p_geometry geometry,
@p_position_as_ratio FLOAT = 0.5,
@p_round_x INT = 3,
@p_round_y INT = 3,
@p_round_z INT = 3
)
RETURNS geometry
AS
BEGIN
DECLARE
@v_geometry geometry;
BEGIN
IF ( @p_geometry IS NULL )
RETURN @p_geometry;
SELECT TOP 1
@v_geometry = geometry::STGeomFromText('POINT(' +
LTRIM(STR(round(i3.x2-((i3.x2-i3.x1)*i3.vectorPositionRatio),@p_round_x),24,15)) + ' ' +
LTRIM(STR(round(i3.y2-((i3.y2-i3.y1)*i3.vectorPositionRatio),@p_round_y),24,15)) +
CASE WHEN i3.z2 IS NOT NULL
THEN ' ' + LTRIM(STR(round(i3.z2-((i3.z2-i3.z1)*vectorPositionRatio),@p_round_z),24,15))
ELSE ''
END + ')',@p_geometry.STSrid)
FROM (SELECT /* select vector/segment "containing" centroid or mid-point of linestring */
i2.SRID,
i2.X1,i2.Y1,i2.Z1,
i2.X2,i2.Y2,i2.Z2,
(i2.cumLength-i2.pointDistance)/i2.vectorLength AS vectorPositionRatio,
CASE WHEN pointDistance BETWEEN CASE WHEN lag(cumLength,1) OVER (ORDER BY VID) IS NULL
THEN 0
ELSE lag(cumLength,1) OVER (ORDER BY VID)
END
AND vectorLength +
CASE WHEN lag(cumLength,1) OVER (ORDER BY VID) IS NULL
THEN 0
ELSE lag(cumLength,1) OVER (ORDER BY VID)
END
THEN 1
ELSE 0
END AS RightSegment
FROM (SELECT i1.VID,
i1.SRID,
i1.X1,i1.Y1,i1.Z1,
i1.X2,i1.Y2,i1.Z2,
i1.vectorLength,
i1.pointDistance,
/* generate cumulative length */
SUM(vectorLength) OVER (ORDER BY vid ROWS UNBOUNDED PRECEDING) AS cumLength
FROM (SELECT v.id AS vid,
@p_geometry.STSrid AS srid,
(@p_geometry.STLength() * @p_position_as_ratio) AS pointDistance,
v.sx AS X1, v.sy AS Y1, v.sz AS Z1,
v.ex AS X2, v.ey AS Y2, v.ez AS Z2,
geometry::Point(v.sX,v.sy/*,v.sz*/,@p_geometry.STSrid).STDistance(geometry::Point(v.eX,v.ey,/*,v.sz*/@p_geometry.STSrid)) AS vectorLength
FROM [dbo].[STVECTORIZE](@p_geometry) AS v
) i1
) i2
) i3
WHERE i3.rightSegment = 1;
RETURN @v_geometry;
END;
END
GO
Testing a range of percentage values (5/10, 0.5 or 50% is centroid)
WITH line AS (
SELECT geometry::STGeomFromText('LINESTRING(258.72254365233152 770.97400259630615, 268.79365642517564 739.08214548229967, 278.86476919801976 707.1902883682933, 332.57737065318844 693.76213800450114, 366.14774656266889 676.97695004976094, 426.57442319973364 697.11917559544918, 520.57147574627891 737.40362668682576, 631.35371624756431 744.11770186872184, 829.41893411349884 797.83030332389046, 1547.8249785763801 791.11622814199438, 1205.4071442996797 895.18439346138371, 832.77597170444687 1039.5370098721496, 490.3581374277465 1086.5355361454222, 416.50331042688953 1076.464423372578, 381.25441572193506 1059.6792354178378, 346.00552101698065 1042.8940474630976, 320.82773908487036 1019.3947843264614, 295.64995715276001 995.89552118982499, 287.25736317538986 964.00366407581862, 278.86476919801976 932.11180696181225, 282.2218067889678 891.82735587043567, 277.18625040254574 858.25697996095528, 272.15069401612368 824.68660405147489, 258.72254365233152 770.97400259630615)',0) AS geom
)
SELECT 'O' AS label, geom FROM line
UNION ALL
SELECT CAST(a.intValue*10 AS VARCHAR) + '%' AS label, [dbo].[STCentroid_L](geom,CAST(a.IntValue AS FLOAT)/10.0,3,3,2).STBuffer(15) AS geom
FROM line
CROSS apply
dbo.generate_series(0,9,1) a;
This looks like:
Generic Centroid function
Selects the right function out of those above.
DROP FUNCTION [dbo].[STCentroid];
CREATE FUNCTION [dbo].[STCentroid] (
@p_geometry geometry,
@p_start INT = 0, /* 0 = use average of all Area's vertices for starting X centroid calculation otherwise centre X of MBR */
@p_largest INT = 1,
@p_round_x INT = 3,
@p_round_y INT = 3,
@p_round_z INT = 3
)
RETURNS geometry
AS
BEGIN
DECLARE
@v_geomn INT,
@v_geom geometry,
@v_val FLOAT,
@v_test_val FLOAT,
@v_geometry geometry,
@v_gtype nvarchar(MAX);
BEGIN
IF ( @p_geometry IS NULL )
RETURN @p_geometry;
SET @v_gtype = @p_geometry.STGeometryType();
SET @v_geometry = @p_geometry;
IF ( @v_gtype IN ('MultiPolygon','MultiLineString') )
BEGIN
SET @v_test_val = -9999999999.9999999;
IF ( @p_largest = 0 )
SET @v_test_val = 9999999999.9999999;
-- Get parts of multi-part geometry
--
SET @v_geomn = 1;
WHILE ( @v_geomn <= @p_geometry.STNumGeometries() )
BEGIN
SET @v_geometry = @p_geometry.STGeometryN(@v_geomn);
IF ( @v_geometry.STGeometryType() = 'LineString' )
SET @v_val = @v_geometry.STLength()
ELSE
SET @v_val = @v_geometry.STArea();
IF ( @p_largest = 1 )
BEGIN
IF ( @v_val >= @v_test_val )
BEGIN
SET @v_test_val = @v_val;
SET @v_geom = @v_geometry;
END;
END;
ELSE
BEGIN
IF ( @v_val <= @v_test_val )
BEGIN
SET @v_test_val = @v_val;
SET @v_geom = @v_geometry;
END;
END;
SET @v_geomn = @v_geomn + 1;
END;
SET @v_geometry = @v_geom;
END;
-- Now get centroid
--
IF ( @v_gtype = 'Point' ) SET @v_geometry = @p_geometry;
IF ( @v_gtype = 'MultiPoint' ) SET @v_geometry = dbo.STCentroid_P(@p_geometry,@p_round_x,@p_round_y,@p_round_z)
IF ( @v_gtype IN ('LineString','MultiLineString' ) ) SET @v_geometry = dbo.STCentroid_L(@v_geometry,0.5,@p_round_x,@p_round_y,@p_round_z)
IF ( @v_gtype IN ('Polygon','MultiPolygon') ) SET @v_geometry = dbo.STCentroid_A(@v_geometry,@p_start,NULL,@p_round_x,@p_round_y,@p_round_z)
RETURN @v_geometry;
END;
END
GO
Testing this we get.
WITH mPoly AS (
SELECT geometry::STGeomFromText('MULTIPOLYGON (((0 0, 100 0, 100 100, 0 100, 0 0)), ((110 110, 150 110, 150 150, 110 150, 110 110)))',0) AS geom
)
SELECT geom FROM mPoly
UNION ALL
SELECT dbo.STCentroid(geom,1 /* Weight By MBR not Vertices */,1 /* Use Largest part in any multipart */,3,3,3).STBuffer(5) AS geom
FROM mPoly;
This looks like:
MultiCentroid function for multi part geometries
DROP FUNCTION [dbo].[STMULTICENTROID];
CREATE FUNCTION [dbo].[STMULTICENTROID] (
@p_geometry GEOMETRY,
@p_start INT = 0, /* 0 = use average of all Area's vertices for starting X centroid calculation otherwise centre X of MBR */
@p_round_x INT = 3,
@p_round_y INT = 3,
@p_round_z INT = 3 )
RETURNS GEOMETRY
AS
BEGIN
DECLARE
@v_geomn INT,
@v_centroid geometry,
@v_part_geom geometry,
@v_geometry geometry,
@v_gtype nvarchar(MAX);
BEGIN
IF ( @p_geometry IS NULL )
RETURN @p_geometry;
SET @v_gtype = @p_geometry.STGeometryType();
IF ( @v_gtype = 'Point' ) SET @v_geometry = @p_geometry;
IF ( @v_gtype = 'LineString') SET @v_geometry = [dbo].[STCENTROID_L](@p_geometry,0.5,@p_round_x,@p_round_y,@p_round_z)
IF ( @v_gtype = 'Polygon' ) SET @v_geometry = [dbo].[STCENTROID_A](@p_geometry,@p_start,NULL,@p_round_x,@p_round_y,@p_round_z)
IF ( @v_gtype = 'MultiPoint' ) SET @v_geometry = @p_geometry;
IF ( @v_gtype IN ('MultiPolygon','MultiLineString') )
BEGIN
-- Get parts of multi-part geometry
--
SET @v_geomn = 1;
WHILE ( @v_geomn <= @p_geometry.STNumGeometries() )
BEGIN
SET @v_part_geom = @p_geometry.STGeometryN(@v_geomn);
IF ( @v_part_geom.STGeometryType() = 'LineString' )
SET @v_centroid = [dbo].[STCENTROID_L](@v_part_geom,0.5,@p_round_x,@p_round_y,@p_round_z)
ELSE
SET @v_centroid = [dbo].[STCENTROID_A](@v_part_geom,@p_start,NULL,@p_round_x,@p_round_y,@p_round_z);
SET @v_geometry = CASE WHEN @v_geometry IS NULL THEN @v_centroid ELSE @v_geometry.STUnion(@v_centroid) END;
SET @v_geomn = @v_geomn + 1;
END;
END;
RETURN @v_geometry;
END;
END
GO
Testing we get.
WITH mPoly AS (
SELECT geometry::STGeomFromText('MULTIPOLYGON (((0 0, 100 0, 100 100, 0 100, 0 0)), ((110 110, 150 110, 150 150, 110 150, 110 110)))',0) AS geom
)
SELECT geom FROM mPoly
UNION ALL
SELECT dbo.STMultiCentroid(geom,1 /* Weight By MBR not Vertices */,3,3,3).STBuffer(5) AS geom
FROM mPoly;
Which looks like:
I hope this is of use to SQL Server Spatial users out there.