Obtención de la longitud máxima del polígono y la anchura media mediante PostGIS

Question

Tengo en PostGIS una tabla tipo polígono.

Necesito calcular automáticamente la longitud máxima del polígono:

Y la segunda es la anchura media:

Estoy dudando si esto es posible ya que, aunque todos mis polígonos serán casi rectangulares en otros casos puede ser realmente ambiguo, y es casi imposible distinguir entre ancho y alto.

Answer 1

Para la primera parte de su pregunta: ¿Qué pasa con ST_LongestLine utilizando la misma geometría dos veces como entrada?

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

Para la segunda parte de su pregunta:

En cuanto al cálculo de la anchura media de los polígonos, aquí se pueden encontrar algunas respuestas interesantes: ¿Cómo puedo calcular la anchura media de un polígono?

Answer 2

Para la primera parte utilice ST_MaxDistance

Devuelve la distancia máxima bidimensional entre dos líneas en unidades proyectadas. Si g1 y g2 es la misma geometría la función devolverá la distancia entre los dos vértices más alejados entre sí en esa geometría.

Exemple :

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

Answer 3

Segunda parte

Para la segunda parte de su pregunta, este viejo gist es una solución: puedes probar la precisión de la solución con tus geometrías "casi rectangulares", y, para "otros casos" te sorprenderá ver que hay una aproximación al rectángulo más cercano. En ambos casos habrá anchura y altura estimaciones.

Uso:

  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

La puesta en práctica:

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;

Esta solución está relacionada con esta discusión .

Primera parte

Para la primera parte de su pregunta, véase La solución de @ThomasB

Answer 4

Así que, si entiendo bien tu pregunta, hay muchas soluciones y ésta es una de ellas.

No queda claro en tu dibujo si quieres tener en cuenta las rotaciones en los objetos poligonales, si todas las formas son similares al patrón, etc.

Para el primer ejemplo (imagen 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)

Para el segundo ejemplo (imagen 2), he creado intuitivamente otra función de cálculo espacial y he calculado en m, la anchura media de un objeto poligonal complejo:

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>

Esta respuesta es la continuación de esa respuesta: https://gis.stackexchange.com/a/418891/120129

Soluciones espaciales originales...

Traducido con www.DeepL.com/Translator (versión gratuita)