¿Calculando el ancho promedio del polígono?

38

Estoy interesado en examinar el ancho promedio de un polígono que representa la superficie de la carretera. También tengo la línea central de la carretera como un vector (que a veces no está exactamente en el centro). En este ejemplo, la línea central del camino está en rojo y el polígono es azul:

Unenfoquedefuerzabrutaenelquehepensadoesamortiguarlalíneaenpequeñosincrementos,intersecarelbúferconunarejilladered,intersectarelpolígonodecarreteraconunarejilladered,calculareláreaintersectadaparaambasmedidasdeintersecciónymantenerhaciendoestohastaqueelerrorseapequeño.Sinembargo,esteesunenfoqueburdo,ymepreguntosihayunasoluciónmáselegante.Además,estoocultaríaelanchodeunacarreteragrandeyunacarreterapequeña.

EstoyinteresadoenunasoluciónqueuseelsoftwareArcGIS10,PostGIS2.0oQGIS.Hevistoesta pregunta y descargué tool para ArcGIS10, pero no pude calcular lo que quiero con él.

Acabo de descubrir la herramienta Geometría de delimitación mínima en ArcGIS 10 que me permite producir los siguientes polígonos verdes:

Esto parece ser una buena solución para carreteras que siguen una cuadrícula, pero no funcionarían de otra manera, por lo que todavía estoy interesado en otras sugerencias.

    
pregunta djq 13.02.2012 - 23:32

4 respuestas

37

Parte del problema es encontrar una definición adecuada de "ancho promedio". Varios son naturales pero diferirán, al menos ligeramente. Para simplificar, considere definiciones basadas en propiedades que son fáciles de calcular (lo que excluirá aquellas basadas en la transformación del eje medial o las secuencias de búferes, por ejemplo).

Como ejemplo, considera que la intuición arquetípica de un polígono con un "ancho" definido es un pequeño búfer (por ejemplo, el radio r con extremos cuadrados) alrededor de un polilínea larga y bastante recta (digamos de longitud L ). Pensamos en 2r = w como su ancho. Así:

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

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

El ancho w y la longitud L se pueden recuperar 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é seguro de que el polígono es realmente largo y delgado, como una aproximación adicional puede tomar 2L + 2w para igualar 2L, desde donde

  • w (crudamente) = 2A / P.

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

Este enfoque no solo es extremadamente simple (solo divide el área por el perímetro y multiplica por 2), con cualquiera de las dos fórmulas, no importa cómo esté orientado el polígono o dónde está situado (porque tales movimientos euclidianos no cambian ni el área ni el perímetro).

Puede considerar usar cualquiera de estas fórmulas para estimar el ancho promedio para cualquier polígono que represente segmentos de calle. El error que 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 curvatura es t radianes (esta es la curvatura absoluta total de la polilínea), entonces realmente

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

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

Conecte estos a la solución anterior (fórmula cuadrática) y simplifique. Cuando se despeja el humo, ¡la contribución del término de curvatura t ha desaparecido! Lo que originalmente parecía una aproximación es perfectamente exacto para los buffers de polilínea que no se intersectan por sí mismos (con extremos cuadrados). Para polígonos de ancho variable, esta es, por lo tanto, una definición razonable de ancho medio.

    
respondido por el whuber 14.02.2012 - 00:47
17

Aquí muestro poca optimización sobre la solución @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 st_buffer, que devuelve una estimación de ancho?

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 calle , sw , la solución es

 sw = buffer_width(ST_Length(g1), g2)

donde sw es el "ancho promedio", g1 la línea central de g2 , y la calle g2 es un POLYGON . Usé solo la biblioteca estándar de OGC, probada con PostGIS , y resolví otras aplicaciones prácticas serias con la misma función de buffer_width.

DEMOSTRACIÓN

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

Suponiendo que podemos generar g2 por g2=ST_Buffer(g1,w) , y que g1 es una recta, entonces g2 es un rectángulo con la longitud L1 y el ancho 2*w , y

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

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

Aquí no evaluamos los búferes con "endcap = square" o "endcap = round", que necesitan una suma para A2 de un área de un búfer de punto con el mismo w .

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

PRUEBAS Y RAZONES

Para las líneas rectas, los resultados, según lo esperado, son exactos. Pero para otras geometrías los resultados pueden ser decepcionantes. La razón principal es que, tal vez, todo el modelo es para rectángulos exactos, o para geometrías que se pueden aproximar a un "rectángulo de franja". Aquí hay un "kit de prueba" para verificar los límites de esta aproximación (vea 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 OTRAS GEOMETRÍAS (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 vea ST_Buffer guide , con buenos ilustratins y los LINESTRING usados aquí.

CONCLUSIONES :

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

Precaución y prevención

¿Por qué se produce el error de estimación? Cuando usa ST_Buffer(g,w) , espera, por el "modelo de tira rectangular", que la nueva área agregada por el búfer de ancho w sea aproximadamente w*ST_Length(g) o w*ST_Perimeter(g) ... Cuando no, generalmente por superposiciones ( ver líneas plegadas) o por "estilo", es cuando la estimación del promedio w fault . Este es el mensaje principal de las pruebas.

Para detectar este problema en cualquier grupo de búferes , compruebe el comportamiento de la generación de búferes:

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%

    
respondido por el Peter Krauss 10.09.2012 - 20:17
13

Si puede unir los datos de sus polígonos a sus datos de la línea central (ya sea por medios espaciales o tabulares), simplemente sume las áreas de los polígonos para cada alineación de la línea central y divídalo por la longitud de la línea central.

    
respondido por el Scro 14.02.2012 - 00:54
8

He desarrollado una fórmula para el ancho promedio de un polígono y la puse en una función Python / ArcPy. Mi fórmula se deriva de (pero se extiende sustancialmente) la noción más directa del ancho promedio que he visto discutida en otra parte; es decir, el diámetro de un círculo que tiene la misma área que tu polígono. Sin embargo, en la pregunta anterior y en mi proyecto, me interesaba más el ancho del eje más estrecho. Además, me interesaba el ancho promedio 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 = ['[email protected]', '[email protected]', 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 el ancho promedio (y algunos otros atributos de geometría para referencia) a través de una variedad de formas usando la función de arriba:

    
respondido por el Tom 23.02.2016 - 01:05

Lea otras preguntas en las etiquetas