# How can I calculate the average width of a polygon?

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:

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:

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.

-
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 – Dan Patterson 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

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, your formula `w = (P - Sqrt(P^2 - 16A)/4` is missing a closing parenthesis. – ForeverWintr Aug 12 '13 at 21:46
@ForeverWintr Thank you for reading the formulas carefully! – whuber Aug 13 '13 at 15:08

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%
``````

-

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

I know this question is old, but in case someone else stumbles across this, 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))
``````

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:

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:

-