¿Calcular la anchura media del polígono?

Question

Me interesa examinar la anchura media de un polígono que representa la superficie de la carretera. También tengo la línea central de la carretera como vector (que a veces no está exactamente en el centro). En este ejemplo, la línea central de la carretera está en rojo y el polígono en azul:

Un enfoque de fuerza bruta que he pensado es amortiguar la línea con pequeños incrementos, intersecar el amortiguador con una rejilla de red, intersecar el polígono de la carretera con una rejilla de red, calcular el área intersecada para ambas medidas de intersección y seguir haciendo esto hasta que el error sea pequeño. Sin embargo, este es un enfoque burdo, y me pregunto si hay una solución más elegante. Además, esto ocultaría la anchura de una carretera grande y otra pequeña.

Estoy interesado en una solución que utilice el software ArcGIS 10, PostGIS 2.0 o QGIS. He visto esto pregunta y descargué el libro de Dan Patterson herramienta para ArcGIS10, pero no he podido calcular lo que quiero con él.

Acabo de descubrir el Geometría de límites mínimos en ArcGIS 10 que me permite producir los siguientes polígonos verdes:

Esto parece una buena solución para las carreteras que siguen una cuadrícula, pero no funcionaría de otro modo, así que sigo interesado en cualquier otra sugerencia.

Answer 1

Parte del problema es encontrar una definición adecuada de "anchura media". Varias son naturales, pero diferirán, al menos ligeramente. Para simplificar, consideremos las definiciones basadas en propiedades fáciles de calcular (lo que va a descartar las basadas en la transformación del eje medio o en las secuencias de topes, por ejemplo).

Como ejemplo, Consideremos que la intuición arquetípica de un polígono con una "anchura" definida es un pequeño tope (digamos de radio r con extremos cuadrados) alrededor de una polilínea larga y bastante recta (digamos de longitud L ). Pensamos en 2r = w como su anchura. Así:

• Su perímetro P es aproximadamente igual a 2L + 2w;

• Su área A es aproximadamente igual a w L.

La anchura w y la longitud L pueden recuperarse entonces como raíces de la cuadrática x^2 - (P/2)x + A; en particular, podemos estimar

• w = (P - Sqrt(P^2 - 16A))/4 .

Cuando estés seguro de que el polígono es realmente largo y delgado, como una aproximación más se puede tomar 2L + 2w para igualar 2L, de donde

• w (groseramente) = 2A / P.

El error relativo de esta aproximación es proporcional a w/L: cuanto más delgado sea el polígono, más se acercará w/L a cero y mejor será la aproximación.

Este enfoque no sólo es extremadamente sencillo (basta con dividir el área por el perímetro y multiplicar por 2), sino que con cualquiera de las dos fórmulas no importa cómo esté orientado el polígono ni dónde esté situado (porque esos movimientos euclidianos no cambian ni el área ni el perímetro).

Puede utilizar cualquiera de estas fórmulas para estimar la anchura media para cualquier polígono que represente segmentos de calle. El error que se comete en la estimación original de w (con la fórmula cuadrática) se produce porque el área A también incluye pequeñas cuñas en cada curva de la polilínea original. Si la suma de los ángulos de las curvas es t radianes (esta es la curvatura total absoluta de la polilínea), entonces realmente

• P = 2L + 2w + 2 Pi t w y

• A = L w + Pi t w^2.

Introduce esto en la solución anterior (fórmula cuadrática) y simplifica. Cuando el humo se disipa, la contribución del término de curvatura t ¡ha desaparecido! Lo que originalmente parecía una aproximación es perfectamente exacto para topes de polilíneas no auto-interseccionados (con extremos cuadrados). Por lo tanto, para los polígonos de anchura variable, se trata de un definición de la anchura media.

Answer 2

Aquí muestro un poco de optimización sobre la solución de @whuber, y estoy poniendo en términos de "ancho de búfer", porque es útil para integrar la solución de un problema más general: ¿Existe una función inversa de st_buffer que devuelva una estimación de la anchura?

CREATE FUNCTION buffer_width(
        -- rectangular strip mean width estimator
    p_len float,   -- len of the central line of g
    p_geom geometry, -- g
    p_btype varchar DEFAULT 'endcap=flat' -- st_buffer() parameter
) RETURNS float AS $f$
  DECLARE
    w_half float;
    w float;    
  BEGIN
         w_half := 0.25*ST_Area(p_geom)/p_len;
         w      := 0.50*ST_Area( ST_Buffer(p_geom,-w_half,p_btype) )/(p_len-2.0*w_half);
     RETURN w_half+w;
  END
$f$ LANGUAGE plpgsql IMMUTABLE;

Para este problema, la pregunta de @celenius sobre ancho de la calle , sw la solución es

 sw = buffer_width(ST_Length(g1), g2)

donde sw es la "anchura media", g1 la línea central de g2 y la calle g2 es un POLÍGONO . He utilizado únicamente la biblioteca estándar OGC, probada con PostGIS , y ha resuelto otras aplicaciones prácticas serias con la misma función buffer_width.

DEMONSTRACIÓN

A2 es el área de g2 , L1 la longitud de la línea central ( g1 ) de g2 .

Suponiendo que podamos generar g2 por g2=ST_Buffer(g1,w) y que g1 es una recta, por lo que g2 es un rectángulo de longitud L1 y la anchura 2*w y

    A2 = L1*(2*w)   -->  w = 0.5*A2/L1

No es la misma fórmula de @whuber, porque aquí w es una mitad del rectángulo ( g2 ) de ancho. Es un buen estimador, pero como podemos ver mediante pruebas (más abajo), no es exacto, y la función lo utiliza como pista, para reducir el g2 y como estimador final.

Aquí no se evalúan los buffers con "endcap=cuadrado" o "endcap=redondo", que necesitan una suma para A2 de un área de un buffer de puntos con el mismo w .

REFERENCIAS: en un foro similar de 2005 W. Huber explica estas y otras soluciones.

PRUEBAS Y RAZONES

Para las líneas rectas los resultados, como era de esperar, son exactos. Pero para otras geometrías los resultados pueden ser decepcionantes. La razón principal es, quizás, que todo el modelo es para rectángulos exactos, o para geometrías que pueden aproximarse a un "rectángulo de tira". Aquí hay un "kit de prueba" para comprobar los límites de esta aproximación (ver wfactor en los resultados anteriores).

 SELECT *, round(100.0*(w_estim-w)/w,1) as estim_perc_error
    FROM (
        SELECT btype, round(len,1) AS len, w, round(w/len,3) AS wfactor,
               round(  buffer_width(len, gbase, btype)  ,2) as w_estim ,
               round(  0.5*ST_Area(gbase)/len       ,2) as w_near
        FROM (
         SELECT
            *, st_length(g) AS len, ST_Buffer(g, w, btype) AS gbase
         FROM (
               -- SELECT ST_GeomFromText('LINESTRING(50 50,150 150)') AS g, -- straight
               SELECT ST_GeomFromText('LINESTRING(50 50,150 150,150 50,250 250)') AS g,
            unnest(array[1.0,10.0,20.0,50.0]) AS w
              ) AS t, 
             (SELECT unnest(array['endcap=flat','endcap=flat join=bevel']) AS btype
             ) AS t2
        ) as t3
    ) as t4;

RESULTADOS:

CON RECTÁNGULOS (la línea central es una LÍNEA RECTA):

         btype          |  len  |  w   | wfactor | w_estim | w_near | estim_perc_error 
------------------------+-------+------+---------+---------+--------+------------------
 endcap=flat            | 141.4 |  1.0 |   0.007 |       1 |      1 |                0
 endcap=flat join=bevel | 141.4 |  1.0 |   0.007 |       1 |      1 |                0
 endcap=flat            | 141.4 | 10.0 |   0.071 |      10 |     10 |                0
 endcap=flat join=bevel | 141.4 | 10.0 |   0.071 |      10 |     10 |                0
 endcap=flat            | 141.4 | 20.0 |   0.141 |      20 |     20 |                0
 endcap=flat join=bevel | 141.4 | 20.0 |   0.141 |      20 |     20 |                0
 endcap=flat            | 141.4 | 50.0 |   0.354 |      50 |     50 |                0
 endcap=flat join=bevel | 141.4 | 50.0 |   0.354 |      50 |     50 |                0

CON OTROS GEOMETROS (línea central doblada):

         btype          | len |  w   | wfactor | w_estim | w_near | estim_perc_error 
 -----------------------+-----+------+---------+---------+--------+------------------
 endcap=flat            | 465 |  1.0 |   0.002 |       1 |      1 |                0
 endcap=flat join=bevel | 465 |  1.0 |   0.002 |       1 |   0.99 |                0
 endcap=flat            | 465 | 10.0 |   0.022 |    9.98 |   9.55 |             -0.2
 endcap=flat join=bevel | 465 | 10.0 |   0.022 |    9.88 |   9.35 |             -1.2
 endcap=flat            | 465 | 20.0 |   0.043 |   19.83 |  18.22 |             -0.9
 endcap=flat join=bevel | 465 | 20.0 |   0.043 |   19.33 |  17.39 |             -3.4
 endcap=flat            | 465 | 50.0 |   0.108 |   46.29 |  40.47 |             -7.4
 endcap=flat join=bevel | 465 | 50.0 |   0.108 |   41.76 |  36.65 |            -16.5

 wfactor= w/len
 w_near = 0.5*area/len
 w_estim is the proposed estimator, the buffer_width function.

Acerca de btype véase Guía de ST_Buffer con buenas ilustraciones y los LINESTRINGs utilizados aquí.

CONCLUSIONES :

• el estimador de w_estim es siempre mejor que w_near ;
• para "cerca de rectangular" g2 geometrías, está bien, cualquier wfactor
• para otras geometrías (cercanas a las "tiras rectangulares"), utilice el límite wfactor=~0.01 para un 1% de error en w_estim . Hasta este factor w, utilice otro estimador.

Precaución y prevención

¿Por qué se produce el error de estimación? Cuando se utiliza ST_Buffer(g,w) se espera, por el "modelo de banda rectangular", que la nueva área añadida por el buffer de ancho w se trata de w*ST_Length(g) o w*ST_Perimeter(g) ... Cuando no es así, normalmente mediante superposiciones (ver líneas dobladas) o mediante "estilización", es cuando la estimación de la media w fallo . Este es el mensaje principal de las pruebas.

Para detectar este problema en cualquier rey del buffer , comprueba el comportamiento de la generación del buffer:

SELECT btype, w, round(100.0*(a1-len1*2.0*w)/a1)::varchar||'%' AS straight_error,  
                 round(100.0*(a2-len2*2.0*w)/a2)::varchar||'%' AS curve2_error,
                 round(100.0*(a3-len3*2.0*w)/a3)::varchar||'%' AS curve3_error
FROM (
 SELECT
    *, st_length(g1) AS len1, ST_Area(ST_Buffer(g1, w, btype)) AS a1,
    st_length(g2) AS len2, ST_Area(ST_Buffer(g2, w, btype)) AS a2,
    st_length(g3) AS len3, ST_Area(ST_Buffer(g3, w, btype)) AS a3
 FROM (
       SELECT ST_GeomFromText('LINESTRING(50 50,150 150)') AS g1, -- straight
              ST_GeomFromText('LINESTRING(50 50,150 150,150 50)') AS g2,
              ST_GeomFromText('LINESTRING(50 50,150 150,150 50,250 250)') AS g3,
              unnest(array[1.0,20.0,50.0]) AS w
      ) AS t, 
     (SELECT unnest(array['endcap=flat','endcap=flat join=bevel']) AS btype
     ) AS t2
) as t3;

RESULTADOS:

         btype          |  w   | straight_error | curve2_error | curve3_error 
------------------------+------+----------------+--------------+--------------
 endcap=flat            |  1.0 | 0%             | -0%          | -0%
 endcap=flat join=bevel |  1.0 | 0%             | -0%          | -1%
 endcap=flat            | 20.0 | 0%             | -5%          | -10%
 endcap=flat join=bevel | 20.0 | 0%             | -9%          | -15%
 endcap=flat            | 50.0 | 0%             | -14%         | -24%
 endcap=flat join=bevel | 50.0 | 0%             | -26%         | -36%

       

Answer 3

Si puede unir los datos de los polígonos con los datos de las líneas centrales (ya sea por medios espaciales o tabulares), sólo tiene que sumar las áreas de los polígonos para cada alineación de la línea central y dividirlas por la longitud de la línea central.

Answer 4

He desarrollado una fórmula para la anchura media de un polígono y la he puesto en una función de Python/ArcPy. Mi fórmula se deriva de (pero amplía sustancialmente) la noción más directa de la anchura media que he visto discutir en otros lugares; es decir, el diámetro de un círculo que tiene la misma área que su polígono. Sin embargo, en la pregunta anterior y en mi proyecto, me interesaba más la anchura del eje más estrecho. Además, me interesaba la anchura media para formas potencialmente complejas y no convexas.

Mi solución fue:

(perimeter / pi) * area / (perimeter**2 / (4*pi))
= 4 * area / perimeter

Eso es:

(Diameter of a circle with the same perimeter as the polygon) * Area / (Area of a circle with the same perimeter as the polygon)

La función es:

def add_average_width(featureClass, averageWidthField='Width'):
    '''
    (str, [str]) -> str

    Calculate the average width of each feature in the feature class. The width
        is reported in units of the feature class' projected coordinate systems'
        linear unit.

    Returns the name of the field that is populated with the feature widths.
    '''
    import arcpy
    from math import pi

    # Add the width field, if necessary
    fns = [i.name.lower() for i in arcpy.ListFields(featureClass)]
    if averageWidthField.lower() not in fns:
        arcpy.AddField_management(featureClass, averageWidthField, 'DOUBLE')

    fnsCur = ['SHAPE@LENGTH', 'SHAPE@AREA', averageWidthField]
    with arcpy.da.UpdateCursor(featureClass, fnsCur) as cur:
        for row in cur:
            perim, area, width = row
            row[-1] = ((perim/pi) * area) / (perim**2 / (4 * pi))
            cur.updateRow(row)

    return averageWidthField

Aquí hay un mapa exportado con la anchura media (y algunos otros atributos de la geometría para la referencia) a través de una variedad de formas utilizando la función de arriba:

Answer 5

Otra solución con eje medial aproximado:

1. Calcular el eje medial aproximado del polígono;
2. Obtenga la longitud del eje medial aproximado;
3. Obtiene la distancia desde ambos extremos del eje hasta el borde del polígono;
4. Suma la longitud del eje y las distancias del paso 3 -- es la longitud aproximada del polígono;
5. Ahora puede dividir el área del polígono por esta longitud y obtener la anchura media del polígono.

El resultado será seguramente erróneo para aquellos polígonos en los que el eje medio aproximado no sea una única línea continua, por lo que se puede comprobar antes del paso 1 y devolver NULL o algo así.

Este es un ejemplo de la función de PostgreSQL (nota: necesita instalar postgis y postgis_sfcgal extensiones):

CREATE FUNCTION ApproximatePolygonLength(geom geometry)
RETURNS float AS $$
    SELECT
        CASE
            /* in case when approximate medial axis is empty or simple line
             * return axis length
             */
            WHEN (ST_GeometryType(axis.axis) = 'ST_LineString' OR ST_IsEmpty(axis.axis))
                THEN axis_length.axis_length
                    + start_point_distance.start_point_distance
                    + end_point_distance.end_point_distance
            /* else geometry is too complex to define length */
            ELSE NULL
        END AS length
    FROM
        LATERAL (
            SELECT
                ST_MakeValid(geom) AS valid_geom
        ) AS valid_geom,
        LATERAL (
            SELECT
                /* `ST_LineMerge` returns:
                 *  - `GEOMETRYCOLLECTION EMPTY`, if `ST_ApproximateMedialAxis` is an empty line (i.e. for square);
                 *  - `LINESTRING ...`, if `ST_ApproximateMedialAxis` is a simple line;
                 *  - `MULTILINESTRING ...`, if `ST_ApproximateMedialAxis` is a complex line
                 *     that can not be merged to simple line. In this case we should return `NULL`.
                 */
                ST_LineMerge(
                    ST_ApproximateMedialAxis(
                        valid_geom.valid_geom
                    )
                ) AS axis
        ) AS axis,
        LATERAL (
            SELECT
                ST_Boundary(valid_geom.valid_geom) AS border
        ) AS border,
        LATERAL (
            SELECT
                ST_Length(axis.axis) AS axis_length
        ) AS axis_length,
        LATERAL (
            SELECT
                ST_IsClosed(axis.axis) AS axis_is_closed
        ) AS axis_is_closed,
        LATERAL (
            SELECT
                CASE WHEN axis_is_closed.axis_is_closed THEN 0
                ELSE
                    ST_Distance(
                        border.border,
                        CASE ST_GeometryType(axis.axis)
                            WHEN 'ST_LineString' THEN ST_StartPoint(axis.axis)
                            /* if approximate medial axis is empty (i.e. for square),
                             * get centroid of geometry
                             */
                            ELSE ST_Centroid(valid_geom.valid_geom)
                        END
                    )
                END AS start_point_distance
        ) AS start_point_distance,
        LATERAL (
            SELECT
                CASE WHEN axis_is_closed.axis_is_closed THEN 0
                ELSE
                    ST_Distance(
                        border.border,
                        CASE ST_GeometryType(axis.axis)
                            WHEN 'ST_LineString' THEN ST_EndPoint(axis.axis)
                            /* if approximate medial axis is empty (i.e. for square),
                             * get centroid of geometry
                             */
                            ELSE ST_Centroid(valid_geom.valid_geom)
                        END
                    )
                END AS end_point_distance
        ) AS end_point_distance;
$$ LANGUAGE SQL;

CREATE FUNCTION ApproximatePolygonWidth(geom geometry)
RETURNS float AS $$
    SELECT
        CASE
            WHEN length IS NULL THEN NULL
            ELSE area.area / length.length
        END AS width
    FROM
        (
            SELECT ApproximatePolygonLength(geom) AS length
        ) AS length,
        (
            SELECT
                ST_Area(
                    ST_MakeValid(geom)
                ) AS area
        ) AS area;
$$ LANGUAGE SQL;

Desventaja:

Esta solución no funcionará con los casos en los que el polígono es casi rectangular y el ser humano puede definir intuitivamente su longitud, pero el eje medio aproximado tiene pequeñas ramas cerca del borde y por lo tanto el algoritmo devuelve Ninguno.

Ejemplo: