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

Question

I have in PostGIS a polygon type table.

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

enter image description here

And the second one is the average width:

enter image description here

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/… — John Powell, 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! — Martin F, 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 — raphael, 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. — Peter Krauss, Jan 20 at 13:27

Fezter · Accepted Answer · 2017-08-21 06:31:50Z

14

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?

edited Aug 21, 2017 at 6:31

Fezter♦

21.2k11 gold badges65 silver badges115 bronze badges

answered Jan 30, 2015 at 17:05

Thomas B

6,3201 gold badge13 silver badges26 bronze badges

Add a comment |

raphael · Accepted Answer · 2015-01-30 17:28:52Z

9

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

answered Jan 30, 2015 at 17:28

raphael

3,20717 silver badges56 bronze badges

Add a comment |

PolyGeo · Accepted Answer · 2022-02-02 23:57:18Z

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

Cyril Mikhalchenko · Accepted Answer · 2021-12-18 15:24:52Z

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)

Stack Exchange Network

current community

your communities

more stack exchange communities

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

4 Answers 4

Second part

First part

Not the answer you're looking for? Browse other questions tagged postgis polygon length or ask your own question.

Linked

Hot Network Questions

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

4 Answers 4

Second part

First part

Not the answer you're looking for? Browse other questions tagged postgis polygon length or ask your own question.

Linked

Related

Hot Network Questions