39

I'm interested in examining the average width of a polygon that represents the road surface. I also have the road centerline as a vector (which is sometimes not exactly in the center). In this example, the road-centerline is in red, and the polygon is blue:

enter image description here

One brute force approach I have thought of, is to buffer the line by small increments, intersect the buffer with a fishnet grid, intersect the road-polygon with a fishnet grid, calculate the intersected area for both intersection measures and to keep doing this until the error is small. This is a crude approach though, and I'm wondering if there is a more elegant solution. In addition, this would conceal the width of a large road and a small road.

I'm interested in a solution that uses ArcGIS 10, PostGIS 2.0 or QGIS software. I have seen this question and downloaded Dan Patterson's tool for ArcGIS10, but I wasn't able to calculate what I want with it.

I've just discovered the Minimum Bounding Geometry tool in ArcGIS 10 which enables me to produce the following green polygons:

enter image description here

This seems like a good solution for roads that follow a grid, but would not work otherwise, so I am still interested in any other suggestions.

closed as too broad by PolyGeo Mar 21 at 8:48

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • Did you rule out possible solutions on the sidebar? ie gis.stackexchange.com/questions/2880/… apparently flagged as a trivial answer to a potentially duplicate post – user681 Feb 13 '12 at 23:05
  • @DanPatterson I didn't see any question like this (many are related, of course). Did you mean that my question was flagged? I didn't understand your second line. – djq Feb 13 '12 at 23:13
  • That related question, @Dan, concerns a different interpretation of "width" (well, actually, the interpretation isn't perfectly clear). The replies there appear to focus on finding width at the widest point, rather than an average width. – whuber Feb 13 '12 at 23:49
  • As @whuber want to centralize here discussions, closing another questions, I suggest that the question is edited to people understand "the estimative of average width of a rectangular strip" – Peter Krauss Sep 10 '12 at 20:29
  • @Peter: As a rectangular strip is a fortiori a polygon, the more general title ought to stand. – whuber Sep 10 '12 at 20:58
38

Part of the problem is to find a suitable definition of "average width." Several are natural but will differ, at least slightly. For simplicity, consider definitions based on properties that are easy to compute (which is going to rule out those based on the medial axis transform or sequences of buffers, for instance).

As an example, consider that the archetypal intuition of a polygon with a definite "width" is a small buffer (say of radius r with squared ends) around a long, fairly straight polyline (say of length L). We think of 2r = w as its width. Thus:

  • Its perimeter P is approximately equal to 2L + 2w;

  • Its area A is approximately equal to w L.

The width w and length L can then be recovered as roots of the quadratic x^2 - (P/2)x + A; in particular, we can estimate

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

When you're sure the polygon really is long and skinny, as a further approximation you can take 2L + 2w to equal 2L, whence

  • w (crudely) = 2A / P.

The relative error in this approximation is proportional to w/L: the skinnier the polygon, the closer w/L is to zero, and the better the approximation gets.

Not only is this approach extremely simple (just divide the area by the perimeter and multiply by 2), with either formula it doesn't matter how the polygon is oriented or where it is situated (because such Euclidean motions change neither the area nor the perimeter).

You might consider using either of these formulas to estimate average width for any polygons that represent street segments. The error you make in the original estimate of w (with the quadratic formula) comes about because the area A also includes tiny wedges at each bend of the original polyline. If the sum of bend angles is t radians (this is the total absolute curvature of the polyline), then really

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

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

Plug these into the previous (quadratic formula) solution and simplify. When the smoke clears, the contribution from the curvature term t has disappeared! What originally looked like an approximation is perfectly accurate for non-self-intersecting polyline buffers (with squared ends). For variable-width polygons, this is therefore a reasonable definition of average width.

  • Thanks @whuber that is a great answer and it has helped me to think about it a lot more clearly. – djq Feb 14 '12 at 0:09
  • @whuber: I am writing a paper and I would need to give a proper ('academic') reference to the method you describe here. Do you have such reference? Does this measure has a name? If not, I can name it after you ! What about the "Huber's width measure" ? – julien Mar 15 '18 at 10:34
  • @julien I don't have any references. This format might work: MISC {20279, TITLE = {Calculating average width of polygon?}, AUTHOR = {whuber (gis.stackexchange.com/users/664/whuber)}, HOWPUBLISHED = {GIS}, NOTE = {URL:gis.stackexchange.com/q/20279/664 (version: 2013-08-13)}, EPRINT = {gis.stackexchange.com/q/20279}, URL = {gis.stackexchange.com/q/20279} } – whuber Mar 15 '18 at 13:43
17

Here I show little optimization about @whuber solution, and I'm putting in terms of "buffer width", because it is useful for integrate the solution of a more general problem: Is there a st_buffer inverse function, that returns a width estimation?

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;

For this problem, the @celenius question about street width, sw, the solution is

 sw = buffer_width(ST_Length(g1), g2)

where sw is the "average width", g1 the central line of g2, and the street g2 is a POLYGON. I used only the OGC standard library, tested with PostGIS, and solved other serious practical applications with the same buffer_width function.

DEMONSTRATION

A2 is the area of g2, L1 the length of the central line (g1) of g2.

Supposing that we can generate g2 by g2=ST_Buffer(g1,w), and that g1 is a straight, so g2 is a rectangle with lenght L1 and width 2*w, and

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

It is not the same formula of @whuber, because here w is a half of the rectangle (g2) width. It is a good estimator, but as we can see by tests (below), is not exact, and the function use it as a clue, to reduce the g2 area, and as a final estimator.

Here we do not evaluate buffers with "endcap=square" or "endcap=round", that need a sum to A2 of an area of a point buffer with the same w.

REFERENCES: in a similar forum of 2005, W. Huber explains suchlike and other solutions.

TESTS AND REASONS

For straight lines the results, as expected, are exact. But for other geometries the results can be disappointing. The main reason is, perhaps, all the model is for exact rectangles, or for geometries that can be approximated to a "strip rectangle". Here a "test kit" for check the limits of this approximation (see wfactor in the results above).

 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;

RESULTS:

WITH RECTANGLES (central line is a STRAIGHT LINE):

         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

WITH OTHER GEOMETRIES (centerline folded):

         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.

About btype see ST_Buffer guide, with good ilustratins and the LINESTRINGs used here.

CONCLUSIONS:

  • the estimator of w_estim is always better than w_near;
  • for "near to rectangular" g2 geometries, its ok, any wfactor
  • for another geometries (near to "rectangular strips"), use the limit wfactor=~0.01 for 1% of error on w_estim. Up to this wfactor, use another estimator.

Caution and prevention

Why the estimation error occurs? When you use ST_Buffer(g,w), you expect, by the "rectangular strip model", that the new area added by the buffer of width w is about w*ST_Length(g) or w*ST_Perimeter(g)... When not, usually by overlays (see folded lines) or by "styling", is when the estimation of average w fault. This is the main message of the tests.

To detect this problem at any king of buffer, check the behavior of the buffer generation:

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;

RESULTS:

         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%

        alert

13

If you can join your polygon data to your centerline data (by either spatial or tabular means), then just sum the polygon areas for each centerline alignment and divide by the length of the centerline.

  • that's true! In this case, my centerlines are not the same length, but I could always join them as one, and split them per polygon. – djq Feb 14 '12 at 0:10
  • If your data is in postgreSQL/postGIS and you have a street id field for centerlines and polygons, then there is no need merge/split, and using aggregate functions your answer is just a query away. I'm slow at SQL, or I would post an example. Let me know if this is how you're are going to solve, and I'll help sort it out (if needed.) – Scro Feb 14 '12 at 0:47
  • Thanks Scro, it's not currently in PostGIS, but it's pretty quick to load in. I think I'll try @whuber's approach first but I will compare it to the results from PostGIS (and thank you for the offer of SQL help, but I should be able to manage). Mainly trying to get the approach clear in my head first. – djq Feb 14 '12 at 2:11
  • +1 This is a nice simple solution for circumstances where it's available. – whuber Feb 14 '12 at 6:01
9

I've developed a formula for the average width of a polygon and put it into a Python/ArcPy function. My formula is derived from (but substantially extends) the most straightforward notion of average width that I've seen discussed elsewhere; that is, the diameter of a circle having the same area as your polygon. However, in the question above and in my project, I was more interested in the width of the narrowest axis. Furthermore, I was interested in the average width for potentially complex, non-convex shapes.

My solution was:

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

That is:

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

The function is:

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

Here is an exported map with the average width (and some other geometry attributes for reference) across a variety of shapes using the function from above:

enter image description here

  • 4
    If you simplify the expression, it will be just area / perimeter * 4. – culebrón Feb 21 '17 at 22:14
  • Thanks, @culebrón. I was going for clarity of concept over simplicity of formula, and I never even thought about simplifying the equation. This should save me some processing time. – Tom Feb 22 '17 at 0:07
0

Another solution with approximate medial axis:

  1. Calculate approximate medial axis of the polygon;
  2. Get length of approximate medial axis;
  3. Get distance from both ends of axis to polygon's border;
  4. Sum axis length and distances from step 3 -- it is approximate polygon's length;
  5. Now you can divide polygon's area by this length and get average width of polygon.

Result will be surely wrong for those polygons where approximate medial axis is not a single continuous line, so you can check it before step 1 and return NULL or something.

examples

Here is example of PostgreSQL function (note: you need to install postgis and postgis_sfcgal extensions):

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;

Disadvantage:

This solution will not work with cases when polygon is almost rectangular and human can intuitevely define its length but approximate medial axis has little branches near the edge and thus algorithm returns None.

Example:

broken example

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