# Getting max length of polygon and average width using PostGIS [closed]

I have in PostGIS a polygon type table.

I need to calculate automatically the maximum length of the polygon:

And the second one is the average width:

I am doubting if this is possible as, although all my polygons will be nearly rectangular in other cases it can be really ambiguous, and it is nearly impossible to distinguish between width and height.

• It's an interesting question, and it definitely is possible. At the very least you simply draw a line from every vertex to every other, polynomial runtime for sure, but a solution. For more elegant solutions, take a look at gis.stackexchange.com/questions/32552/… Jan 30, 2015 at 16:24
• Your illustration of "max length" seems clear, but your example "avg width" is not perpendicular to said length. Shouldn't it be? Then you mention "height". All a bit confusing! Jan 30, 2015 at 16:30
• I think more detail is needed on this concept of average width, in particular with examples of the actual geometries you are using Jan 30, 2015 at 17:30
• Why not best fit rectangle and its diagonal? (instead the maximum length) ... Imagine a rectangle with same area and a width parallel to the ST_OrientedEnvelope or similar thing. So, the "average and best fit rectangle", with width W and length L, will supply a diagonal D, `D=SQRT(W^2+D^2)`.... If it is valid, you can use my answer below for all values, W, L and D. Jan 20 at 13:27

For the first part of your question: What about ST_LongestLine using the same geometry twice as input?

``````SELECT
ST_Length(ST_LongestLine(
(SELECT geom FROM mylayer WHERE gid=1),
(SELECT geom FROM mylayer WHERE gid=1))
);
``````

For the second part of your question:

Concerning the calculation of the average width of polygons some interesting answers can be found here: How can I calculate the average width of a polygon?

For part one use ST_MaxDistance

Returns the 2-dimensional maximum distance between two linestrings in projected units. If g1 and g2 is the same geometry the function will return the distance between the two vertices most far from each other in that geometry.

Example:

``````SELECT
gid,
ST_MaxDistance(geom, geom) AS "Max Length"
FROM layer
``````

## Second part

For the second part of your question, this old gist is a solution: you can test the precision of the solution with your "nearly rectangular" geometries, and, for "other cases" you will be surprised to see that there is an approximation to the closest rectangle. In both cases there will be width and height estimations.

Usage:

``````  SELECT shapedescr_sizes(geom) as array_descriptor;
-- Average width is array_descriptor[1]
-- Average height is array_descriptor[2]

-- or
SELECT shapedescr_sizes_tr(shapedescr_sizes(geom)) as report;
-- Human-readable report for array_descriptor
``````

Implementation:

``````CREATE FUNCTION shapedescr_sizes(
-- Shape-descriptor "as rectangle" for geometry description by sizes.
gbase geometry,              -- input
-- p_seqs integer DEFAULT 8,    -- constant constant
-- p_shape varchar DEFAULT '',  -- endcap indicator
p_decplacesof_zero integer DEFAULT 6, -- precision of zero when rounding delta
p_dwmin float DEFAULT 99999999.0,     -- change to ex. 0.0001, if to use.
p_deltaimpact float DEFAULT 9999.0    -- internal constant
)  RETURNS float[] AS \$f\$
DECLARE
ret float[];
dw float;
b float;
L_estim float; -- L as width or length
H_estim float; -- H as height
aorig float;
gaux geometry;
g1 geometry;
A0 float;
A1 float;
c float;
delta float;
per float;
errcod float;
BEGIN
errcod=0.0;
IF gbase IS NULL OR NOT(ST_IsClosed(gbase)) THEN
errcod=1;                  -- ERROR1 (die)
RAISE EXCEPTION 'error %: invalid input geometry',errcod;
END IF;
A0 := ST_Area(gbase);
per := st_perimeter(gbase);
dw := sqrt(A0)/p_deltaimpact;
IF dw>p_dwmin THEN dw:=p_dwmin; END IF;
g1 = ST_Buffer(gbase,dw);
A1 = ST_area(g1);
IF A0>A1 THEN
errcod=10;                 -- ERROR2 (die)
RAISE EXCEPTION 'error %: invalid buffer/geometry with A0=% g.t. A1=%',errcod,A0,A1;
END IF;
IF (A1-A0)>1.001*dw*per THEN
gaux := ST_Buffer(g1,-dw);  -- closing operation.
A0 = ST_Area(gaux);         -- changed area
per := ST_Perimeter(gaux);  -- changed
errcod:=errcod + 0.1;       -- Warning3
END IF;
C := 2.0*dw;
b := -(A1-A0)/C+C;
delta := b^2-4.0*A0;
IF delta<0.0 AND round(delta,p_decplacesof_zero)<=0.0 THEN
delta=0.0; -- for regular shapes like the square
errcod:=errcod + 0.01;  -- Warning2
END IF;
IF delta<0.0 THEN
L_estim := NULL;
H_estim := NULL;
errcod:=errcod+100;        -- ERROR3
ELSE
L_estim := (-b + sqrt(delta))/2.0;
H_estim := (-b - sqrt(delta))/2.0;
END IF;
IF abs(A0-L_estim*H_estim)>0.001 THEN
errcod:=errcod + 0.001;    -- Warning1
END IF;
ret := array[L_estim,H_estim,a0,per,dw,errcod];
return ret;
END
\$f\$ LANGUAGE plpgsql IMMUTABLE;

CREATE or replace FUNCTION shapedescr_sizes_tr(
-- Human translator for shapedescr_sizes(). Uses ROUND(float).
float[],               -- shapedescr_sizes() returned vector
integer DEFAULT 0,     -- general round parameter
integer DEFAULT 3      -- number of "decimal warnings"
) RETURNS varchar[] -- length, height, area, perimeter, dw, radius, err_message
AS \$f\$
SELECT array[
round(L,\$2+1)::varchar, round(H,\$2+1)::varchar, round(area,\$2)::varchar,
round(perim,\$2+1)::varchar, round(dw,\$2+3)::varchar,
round(sqrt(L*H/pi()),\$2+1)::varchar,  -- radius for "shape as circle"
CASE WHEN errcod>3.0 THEN 'ERROR '  ||round(errcod-\$3)
WHEN errcod>0.0  THEN 'WARNING '||round(errcod)||CASE
WHEN round(10^(-errcod),\$3-1)!=(\$1)[6] THEN '.'|| (\$1)[6]*10^\$3
ELSE ''
END
ELSE ''
END::varchar
]
FROM (
SELECT (\$1)[1] as L, (\$1)[2] as H, (\$1)[3] as area, (\$1)[4] as perim,
(\$1)[5] as dw, log((\$1)[6]*10.0^(\$3+1)+1.0) as errcod
) as t;
\$f\$ LANGUAGE SQL IMMUTABLE;
``````

This solution is related with this discussion.

## First part

For the first part of your question see @ThomasB's solution

So, if I understand your question correctly, there are many solutions and this is one of them.

It is not clear from your drawing whether you want to take into account rotations in polygonal objects, whether all shapes are similar to the pattern, etc.

For the first example (picture 1), RUN CTE

``````WITH
tbla AS (SELECT (ST_Dump(geom)).geom geom FROM <table_name>),
tblb AS (SELECT ST_LongestLine(a.geom, b.geom) geom FROM tbla a JOIN
tbla b ON ST_Intersects(a.geom, b.geom))
(SELECT geom, ST_Length(geography(geom)) lenght FROM tblb)
``````

For the second example (picture 2), I intuitively created another spatial calculation function and calculated in m, the average width of a complex polygonal object:

``````CREATE OR REPLACE FUNCTION ST_AvgLengthLineIRRegularPolygon(
geom GEOMETRY,
n integer)
RETURNS double precision AS
\$BODY\$
WITH
tbl_rigth AS (WITH
tbla AS (SELECT (ST_Dump(geom)).geom geom),
tblb AS (SELECT ST_LongestLine(a.geom, b.geom) geom FROM tbla a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
tblc AS (SELECT geom, ST_Length(geography(geom)) lenght FROM tblb),
tbld AS (SELECT ST_Buffer(geom, lenght/1855/60) geom FROM tblc),
tble AS (SELECT ST_Boundary(ST_OrientedEnvelope(geom)) geom FROM tbld),
tblf AS (SELECT line1, line2 FROM (SELECT ST_MakeLine(ST_PointN(geom,1), ST_PointN(geom,2)) line1,
ST_MakeLine(ST_PointN(geom,4), ST_PointN(geom,3)) line2 FROM tble) foo),
tblg AS (SELECT generate_series (0, n) as steps),
tblh AS (SELECT steps AS stp1, ST_LineInterpolatePoint(line1, steps/(SELECT count(steps)::float-1 FROM tblg)) geom1 FROM tblf, tblg GROUP BY tblg.steps, geom1),
tbli AS (SELECT steps AS stp2, ST_LineInterpolatePoint(line2, steps/(SELECT count(steps)::float-1 FROM tblg)) geom2 FROM tblf, tblg GROUP BY tblg.steps, geom2),
tblj AS (SELECT ST_MakeLine(geom1, geom2) geom FROM tblh JOIN tbli ON true AND stp1=stp2),
tblk AS (SELECT (ST_Dump(ST_Intersection(a.geom, b.geom))).geom geom, count(*) cnt FROM tblj a JOIN tbla b ON true GROUP BY a.geom, b.geom)
(SELECT ST_Length(geography(ST_Union(geom)))/SUM(cnt) avg_length FROM tblk)),
tbl_lefth AS (WITH
tbla AS (SELECT (ST_Dump(geom)).geom geom),
tblb AS (SELECT ST_LongestLine(a.geom, b.geom) geom FROM tbla a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
tblc AS (SELECT geom, ST_Length(geography(geom)) lenght FROM tblb),
tbld AS (SELECT ST_Buffer(geom, lenght/1855/60) geom FROM tblc),
tble AS (SELECT ST_Boundary(ST_OrientedEnvelope(geom)) geom FROM tbld),
tblf AS (SELECT line1, line2 FROM (SELECT ST_MakeLine(ST_PointN(geom,2), ST_PointN(geom,3)) line1,
ST_MakeLine(ST_PointN(geom,1), ST_PointN(geom,4)) line2 FROM tble) foo),
tblg AS (SELECT generate_series (0, n) as steps),
tblh AS (SELECT steps AS stp1, ST_LineInterpolatePoint(line1, steps/(SELECT count(steps)::float-1 FROM tblg)) geom1 FROM tblf, tblg GROUP BY tblg.steps, geom1),
tbli AS (SELECT steps AS stp2, ST_LineInterpolatePoint(line2, steps/(SELECT count(steps)::float-1 FROM tblg)) geom2 FROM tblf, tblg GROUP BY tblg.steps, geom2),
tblj AS (SELECT ST_MakeLine(geom1, geom2) geom FROM tblh JOIN tbli ON true AND stp1=stp2),
tblk AS (SELECT (ST_Dump(ST_Intersection(a.geom, b.geom))).geom geom, count(*) cnt FROM tblj a JOIN tbla b ON true GROUP BY a.geom, b.geom)
(SELECT ST_Length(geography(ST_Union(geom)))/SUM(cnt) avg_length FROM tblk))
SELECT avg_length FROM tbl_rigth UNION SELECT avg_length FROM tbl_lefth
\$BODY\$
LANGUAGE SQL
``````

RUN

``````SELECT ST_AvgLengthLineIRRegularPolygon(geom, 300) avglengt FROM <table_name>
``````

This answer is a continuation of that answer: https://gis.stackexchange.com/a/418891/120129

Original spatial solutions...

Translated with www.DeepL.com/Translator (free version)