Je suis intéressé à examiner la largeur moyenne d'un polygone qui représente la surface de la route. J'ai également la ligne centrale de la route comme vecteur (qui parfois n'est pas exactement au centre). Dans cet exemple, l'axe de la route est en rouge et le polygone est en bleu:

Une approche de force brute à laquelle j'ai pensé consiste à tamponner la ligne par petits incréments, à intersecter le tampon avec un maillage en filet, à intersecter le polygone routier avec un maillage à mailles, à calculer la surface intersectée pour les deux mesures d'intersection et à continuer ainsi l'erreur est petite. C'est une approche rudimentaire cependant, et je me demande s'il existe une solution plus élégante. En outre, cela dissimulerait la largeur d’une grande route et d’une petite route.

Je suis intéressé par une solution utilisant les logiciels ArcGIS 10, PostGIS 2.0 ou QGIS. J'ai vu cette question et téléchargé l' outil de Dan Patterson pour ArcGIS10, mais je n'ai pas pu calculer ce que je veux avec.

Je viens de découvrir l' outil de géométrie de limitation minimale dans ArcGIS 10 qui me permet de produire les polygones verts suivants:

Cela semble être une bonne solution pour les routes qui suivent une grille, mais ne fonctionnerait pas autrement, alors je suis toujours intéressé par toute autre suggestion.

Une partie du problème consiste à trouver une définition appropriée de «largeur moyenne». Plusieurs sont naturels mais différeront, au moins légèrement. Pour plus de simplicité, considérons des définitions basées sur des propriétés faciles à calculer (ce qui exclut celles basées sur la transformation de l’axe médian ou sur des séquences de tampons, par exemple).

A titre d'exemple, considérons que l'intuition archétypale d'un polygone ayant une "largeur" ​​définie est un petit tampon (disons de rayon r aux extrémités carrées) autour d'une polyligne longue et assez droite (disons de longueur L ). Nous pensons que 2r = w est sa largeur. Ainsi:

• Son périmètre P est approximativement égal à 2L + 2w;

• Sa surface A est approximativement égale à w L.

La largeur w et la longueur L peuvent alors être récupérées en tant que racines du quadratique x ^ 2 - (P / 2) x + A; en particulier, on peut estimer

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

Lorsque vous êtes certain que le polygone est vraiment long et maigre, vous pouvez prendre une autre approximation de 2L + 2w pour obtenir 2L.

• w (grossièrement) = 2A / P.

L'erreur relative dans cette approximation est proportionnelle à w / L: plus le polygone est mince, plus w / L est proche de zéro et meilleure est l'approximation.

Non seulement cette approche est-elle extrêmement simple (il suffit de diviser l'aire par le périmètre et de la multiplier par 2), quelle que soit la formule choisie, l'orientation du polygone ou son emplacement ne sont pas gênants (car de tels mouvements euclidiens ne changent ni l'aire ni le périmètre).

Vous pouvez envisager d'utiliser l'une de ces formules pour estimer la largeur moyenne de tout polygone représentant des segments de rue. L'erreur que vous faites dans l'estimation initiale de w (avec la formule quadratique) provient du fait que la zone A comprend également de minuscules coins à chaque pli de la polyligne originale. Si la somme des angles de courbure est t radians (c’est la courbure totale absolue de la polyligne), alors

• P = 2L + 2w + 2 Pi tw et

• A = L w + Pi tw ^ 2.

Branchez-les dans la solution précédente (formule quadratique) et simplifiez-les. Lorsque la fumée se dissipe, la contribution du terme de courbure t a disparu! Ce qui ressemblait à l’origine à une approximation est parfaitement exact pour les tampons polylignes ne se croisant pas automatiquement (avec des extrémités carrées). Pour les polygones de largeur variable, il s'agit donc d'une définition raisonnable de la largeur moyenne.

Ici, je montre peu d'optimisation à propos de la solution @whuber, et je parle de "largeur de tampon", car cela est utile pour intégrer la solution d'un problème plus général: existe-t-il une fonction inverse de st_buffer, qui renvoie une estimation de la largeur?

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;

Pour ce problème, la question @celenius au sujet de la largeur de la rue , swla solution est

 sw = buffer_width(ST_Length(g1), g2)

où swest la "largeur moyenne", g1la ligne centrale de g2, et la rue g2est un POLYGONE . J'ai utilisé uniquement la bibliothèque standard OGC, testée avec PostGIS , et résolu d'autres applications pratiques sérieuses avec la même fonction buffer_width.

MANIFESTATION

A2est l'aire de g2, L1la longueur de la ligne centrale ( g1) de g2.

Supposons que nous puissions générer g2par g2=ST_Buffer(g1,w), et qu’il g1s’agisse d’une droite, ainsi que d’ g2un rectangle de longueur L1et de largeur 2*w, et

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

Ce n'est pas la même formule de @whuber, car voici wla moitié de la g2largeur de rectangle ( ). C'est un bon estimateur, mais comme nous pouvons le constater par les tests (ci-dessous), ce n'est pas exact, et la fonction l'utilise comme indice, pour réduire la g2surface, et comme estimateur final.

Ici, nous n'évaluons pas les tampons avec "endcap = square" ou "endcap = round", qui nécessitent la somme A2 d'une surface d'une mémoire tampon de points avec le même w.

RÉFÉRENCES: dans un forum similaire de 2005 , W. Huber explique des solutions similaires .

TESTS ET RAISONS

Pour les lignes droites, les résultats, comme prévu, sont exacts. Mais pour d'autres géométries, les résultats peuvent être décevants. La raison principale en est que, peut-être, tout le modèle est-il pour des rectangles exacts, ou pour des géométries pouvant être approchées à un "rectangle en bande". Voici un "kit de test" pour vérifier les limites de cette approximation (voir wfactorles résultats ci-dessus).

 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;

RÉSULTATS:

AVEC RECTANGLES (la ligne centrale est une ligne droite):

         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

AVEC AUTRES GEOMETRIES (ligne centrale pliée):

         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.

À propos, btypevoir le guide ST_Buffer , avec les bons ilustratins et les LINESTRING utilisés ici.

CONCLUSIONS :

• l'estimateur de w_estimest toujours meilleur que w_near;
• pour les g2géométries "presque rectangulaires" , c'est ok, toutwfactor
• pour les autres géométries (proches des "bandes rectangulaires"), utilisez la limite wfactor=~0.01de 1% d'erreur sur w_estim. Jusqu'à ce facteur, utilisez un autre estimateur.

Attention et prévention

Pourquoi l'erreur d'estimation se produit? Lorsque vous utilisez ST_Buffer(g,w), prévu, par le « modèle de bande rectangulaire », que la nouvelle zone ajoutée par le tampon de la largeur west d' environ w*ST_Length(g)ou w*ST_Perimeter(g)... Lorsqu'ils ne sont pas, en général par des superpositions (voir les lignes pliées) ou en « coiffant », est quand l'estimation de la wfaute moyenne . C'est le message principal des tests.

Pour détecter ce problème à n’importe quel roi de tampon , vérifiez le comportement de la génération de tampon:

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;

RÉSULTATS:

         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%

Si vous pouvez joindre vos données de polygone à vos données de trait d'axe (par des moyens spatiaux ou tabulaires), additionnez simplement les zones de polygone pour chaque alignement de trait d'axe et divisez-les par la longueur du trait d'axe.

J'ai développé une formule pour la largeur moyenne d'un polygone et je l'ai insérée dans une fonction Python / ArcPy. Ma formule est dérivée de (mais étend substantiellement) la notion la plus simple de largeur moyenne dont j'ai déjà discuté ailleurs; c'est-à-dire le diamètre d'un cercle ayant la même surface que votre polygone. Cependant, dans la question ci-dessus et dans mon projet, je m'intéressais davantage à la largeur de l'axe le plus étroit. De plus, je m'intéressais à la largeur moyenne des formes potentiellement complexes non convexes.

Ma solution était:

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

C'est:

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

La fonction est:

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

Voici une carte exportée avec la largeur moyenne (et quelques autres attributs de géométrie à titre de référence) pour diverses formes à l'aide de la fonction précédente:

Une autre solution avec axe médian approximatif:

1. Calculez l’axe médial approximatif du polygone;
2. Obtenir la longueur de l'axe médial approximatif;
3. Obtenez la distance entre les deux extrémités de l'axe et la frontière du polygone;
4. La somme de la longueur de l'axe et des distances depuis l'étape 3 est la longueur approximative du polygone;
5. Vous pouvez maintenant diviser la surface du polygone par cette longueur et obtenir la largeur moyenne du polygone.

Le résultat sera sûrement faux pour les polygones où l’axe médian approximatif n’est pas une simple ligne continue. Vous pouvez donc le vérifier avant l’étape 1 et renvoyer NULLun résultat.

Voici un exemple de fonction PostgreSQL (note: vous devez installer les extensions postgis et postgis_sfcgal ):

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;

Désavantage:

Cette solution ne fonctionnera pas dans les cas où le polygone est presque rectangulaire et que l’être humain peut définir sa longueur de manière exhaustive, mais que l’axe médial approximatif a de petites branches proches du bord et que l’algorithme ne renvoie aucun.

Exemple: