# Is there a st_buffer inverse function, that returns a width estimation?

When you supply a geometry `g1` and a width `w` to the standard OGC function ST_Buffer, it returns a new geometry `g2`,

``````  g2 = ST_Buffer(g1,w)
``````

This question is about the inverse: a function where you can supply two related geometries — `g1` and `g2` satisfies ST_Within(g1,g2) — and returns the width parameter,

``````  w = Unbuffer(g1,g2)
``````

Now the interpretation of `w` is "average width". So, where is this function Unbuffer defined? Are there references for it? Is there a standard library for it?

Complementing after bounty.

Examples of applications (estimation of): street width, sidewalk width, river width, canal width, internal polygon semelhance, riparian vegetation buffer width, etc.

Example of proposals of credible sources: Pavelsky & Smith 2008 ; ArcMap Script Procedure. They are not solutions here because are (very) CPU consumer, "river-specialized", and can not be implemented by OGC standards. They calculate average of river widthes by each of many perpendicular (of a central line) lines; so, they are "brute force" algorithms.

Example of proposal without credible and/or official sources: troubleshooter here. In this case, you can get the bounty finding sources that makes credible the posted solution.

• A possible answer, if no other post a better one here, is here. – Peter Krauss Sep 12 '12 at 10:03

I'm not an expert in Postgis, but I thinks this might work:

``````ST_Distance(ST_Boundary(g2),g1)
``````

The ST_Boundary will transform your Polygon into a Linestring, If g2 was the result of a X distance buffer from g1, then all the vertices of the resulting linestring must be at the same distance X of g1, and ST_Distance (minimun distance) will give you the desired result.

NOTE: This won't work with buffers generated from ST_Buffer with join=benvel or encap=flat.

To resolve the problem with the benvel and the flat. Here my two cents, you could transform the g1 linestring or polygon to the n line segments (Ls) that compose it with this answer, them calculate for each of them the distance of its middle points to g2 with:

``````ST_Distance(ST_Boundary(g2),line_interpolate_point(Ls,0.5))
``````

And finally calculate the Mode of all the those distances.

Assuming that your g1 geometry is not too complex, the mode will give you the exact width of the original buffer, while the mean will give always an aproximation.

• I appreciate your geometric reasoning, your `ST_Distance(ST_Boundary(g2),g1)` solution is useful as a reliable measure of minimum width. But the question is about "average width"... Your idea can be translated/optimized as the `(ST_perimeter(g2) - ST_perimeter(g1))/8.0` estimator, when g1 and g2 are rectangular. A complex calculus, using average of minimum distantes was used here. – Peter Krauss Sep 15 '12 at 23:21
• For a non benvel ou flat encap buffer the minimum width will be the same as the exact width, and the same as the average width. I'm not sure if I understand what is your final objective for this question, clarifying that might help to improve the answers. Is the g2 geometry really a result of a buffer operation? Or is an independent polygon from g1 and you want to find out the most aproximate buffer of g1, by finding the w average width? What do you mean by deal with endcap/benvel join? Can you show\draw the worst case scenario of g1 and g2 geometries? – Alexandre Neto Sep 16 '12 at 0:41
• About use of split line, ok, please complete your answer with a complete function, is possible as a "brute force algorithm" (!). Please also post sources/references, if you have... But, before, check if I understand: if you simulate with ST_Buffer, and test with your minimum distance, of course, it is the same value than average. But check real cases, see the river width application example (see my refs) or the real case of street width, or any other... The real life need averages (!). – Peter Krauss Sep 16 '12 at 14:59

you can try to get vertices from the polygon (buffer result) , using the vertices as input for voronoi polygons , these polygons draws the mid line between the vertices, then try to extract the mid line, by erasing the lines which touches the (buffer result) contour.

i wish it could help

• I understand that for PostGIS the Voronoi function is only a item of WishList... trac.osgeo.org/postgis/ticket/1215 – Peter Krauss Sep 12 '12 at 21:01
• so sorry, i've never used tried this in postgis, but this link could help you to run voronoi polygons in postgis geogeek.garnix.org/2010/09/voronoi-diagrams-in-postgis.html – geogeek Sep 12 '12 at 21:08
• Ops, sorry truncated comments. Thanks about links, I need do some homework with them. – Peter Krauss Sep 12 '12 at 21:11
• search for this paragraph "Thiessen polygon coverage" in this link portal.opengeospatial.org/modules/admin/… – geogeek Sep 12 '12 at 21:16
• Hum.... Ok, now I understand a bit more of your answer... Perhaps you have confused the "geometric inverse" with "inverse for obtain the radius". In the question we not need to generate a central line geometry (we have it as g1 parameter), or skeleton, we need only a float number that represents the "average width", that was the buffer radius (that produced an approximation to g2). – Peter Krauss Sep 12 '12 at 21:53

After ask, discuss, pay bounty, and post other variants of the same question... I decided to build my solution. It is not so generic as I want, but I am working... There are 3 parts, corresponding to the 3 alternative kinds of `g1` inputs: point, line or polygon.

# Solution for point-buffer

From my question, about the width of the beach of a little island polygon.

``````CREATE FUNCTION pointBuffer_width(g geometry, N integer DEFAULT 8) RETURNS float AS \$\$
-- Returns the value of W in the approximation "g ~ ST_Buffer(ST_Centroid(g),W,'quad_segs=N')"
-- where N>0 indicates the reference-shape, 1=square, 2=heptagon, ... circle.
--
SELECT sqrt(  area/(0.5*n*sin(2*pi()/n))  )
FROM ( SELECT ST_Area(\$1) as area,
(CASE WHEN \$2<=1 THEN 4 ELSE 4*\$2 END)::float as n
) as t;
\$\$ LANGUAGE SQL immutable;
``````

# Solution for line-buffer

From @celenius question, about the average width of a street.

``````CREATE FUNCTION lineBuffer_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;
``````

# Solution for polygon-buffer

From my first question about buffers, about estimation of the average width of a sidewalk...

``````CREATE FUNCTION polyBuffer_width(geometry, geometry) RETURNS float AS \$\$
SELECT 2.0*(ST_Area(\$2)-ST_Area(\$1))/(st_perimeter(\$2)+st_perimeter(\$1));
\$\$ LANGUAGE SQL immutable;
``````

# Putting it all together

(... WORKING ...)

``````CREATE or replace FUNCTION buffer_width(
-- Mean width estimator. Returns w_estim for approximates gbase~ST_Buffer(gref,w_estim)
gbase geometry, -- polygon that ST_Contains(gbase,gref)
gref geometry DEFAULT NULL, -- point or line or NULL.
len1 float DEFAULT 0.0, -- len gref, 0 if point, NULL if polygon.
N integer DEFAULT 8,  -- buffer parameter, quad_segs=N.
bufferType varchar DEFAULT 'endcap=flat' -- st_buffer() parameter
) RETURNS float AS \$f\$
DECLARE
w_half float;
w float;
aux float;
greftype varchar;
BEGIN
greftype := geometrytype(gbase);
IF greftype LIKE '%LINE%' OR len1>0.0 THEN
IF gref is not NULL THEN
len1 = ST_Length(gref);
END IF;
w_half := 0.25*ST_Area(gbase)/len1;
w := 0.50*ST_Area( ST_Buffer(gbase,-w_half,bufferType) )/(len1-2.0*w_half);
RETURN w_half+w;
ELSEIF greftype LIKE '%POINT%' OR gref IS NULL THEN
aux = (CASE WHEN N<=1 THEN 4 ELSE 4*N END)::float;
RETURN sqrt(  ST_Area(gbase)/(0.5*aux*sin(2.0*pi()/aux))  );
ELSEIF greftype LIKE '%POLY%' THEN
RETURN 2.0*(ST_Area(gbase)-ST_Area(gref))/(st_perimeter(gbase)+st_perimeter(gref));
ELSE
RETURN -1.0; -- ERROR FLAG.
END IF;
END
\$f\$ LANGUAGE plpgsql IMMUTABLE;
``````

``````CREATE FUNCTION buffer_width(  -- façade for point
gbase geometry,
N integer
) RETURNS float AS \$f\$
SELECT buffer_width(\$1,NULL,0.0,\$2,'');
\$f\$ LANGUAGE SQL IMMUTABLE;

CREATE FUNCTION buffer_width(  -- façade for line
len1 float,
gbase geometry,
bufferType varchar DEFAULT 'endcap=flat'
) RETURNS float AS \$f\$
SELECT buffer_width(\$2,NULL,\$1,8,\$3);
\$f\$ LANGUAGE SQL IMMUTABLE;
``````

# Testing

``````-- TEST-KIT FOR BUFFER-POINT
SELECT w,
(array['square','octagon','','','','','','~circle'])[shape] as shape,
round(buffer_width(circle8,shape),1) as circle8,
round(buffer_width(octagon,shape),1) as octagon,
round(buffer_width(square,shape),1) as square
FROM (
SELECT
w, shape, -- 1, 2, 3 or 10
ST_Buffer(g, w) as circle8, -- default 8, use 20 for perfect circle
FROM ( SELECT ST_GeomFromText('POINT(100 90)') as g,
unnest(array[1.0,150.0]) as w,
unnest(array[1,2,8]) as shape
) as t
) as t2;

-- TEST-KIT FOR BUFFER-LINE
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;
``````