Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
    
A possible answer, if no other post a better one here, is here. –  Peter Krauss Sep 12 '12 at 10:03
add comment

3 Answers 3

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.

share|improve this answer
    
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
add comment

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

share|improve this answer
    
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
1  
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
show 2 more comments
up vote 0 down vote accepted

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; 

Overloads:

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
    ST_Buffer(g, w, 'quad_segs=2') as octagon,
    ST_Buffer(g, w, 'quad_segs=1') as square
 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;
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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