Whenever possible I recently try to get all my GIS work done in QGIS. Most of the time this is no problem at all. Sometimes it makes things even easier, such as when you’re trying to work with your geospatial data in a PostgreSQL/PostGIS database (good luck trying that in ArcGIS!). But sometimes you come across a task that is just so exotic that nobody has ever come across it. Or at least nobody wrote about coming across it…
In my current case I suppose the latter happened, since I can’t imagine being the first QGIS user ever to have the need for concentric ring buffers (a.k.a “donut buffers”). As you might have heard the new QGIS 2.4 came out this weekend, but unfortunately I didn’t have enough time to install it (let alone both of my machines being currently employed in productive work, leaving no backup machine for trying out fresh-from-the-compiler software…) so I had to use the “old” 2.0 version. There might as well be some added functionality in this regard in the latest version, but to my best knowledge QGIS 2.0 does not allow for the generation of multiple buffers around the same feature(s) in one go, let alone multiple ring buffers.
But fear not, that’s exactly one of the reasons why open source software is so great: you can just extend it to fit your needs. In this specific case I opted for a solution inside the database, since both the input and output data would live there anyways. So I was looking for a solution in PostgreSQL/PostGIS.
Inspired by a great blog post by Egge-Jan Pollé I then came up with the following solution in pure SQL and PL/pgSQL:
First we need to prepare the output table
DROP TABLE IF EXISTS schema.tablename; CREATE TABLE schema.tablename ( gid serial PRIMARY KEY, description character varying(50), radius integer, geom geometry(POLYGON, 4326) ); CREATE INDEX idx_geom_schema_tablename ON schema.tablename USING GIST(geom);
Obviously, the schema, table name, field names and types as well as the index name can be chosen at your liking. The same goes for the SRS, in my case the data is using WGS84 (EPSG 4326).
Then follows the interesting part, the PL/pgSQL function that creates the actual buffers:
DROP FUNCTION IF EXISTS concentric_ring_buffers(minimum integer, maximum integer, step integer); CREATE FUNCTION concentric_ring_buffers(minimum integer, maximum integer, step integer) RETURNS integer AS $$ DECLARE quantity integer := 0; BEGIN RAISE NOTICE 'Generating concentric ring buffers from % km to % km in % km steps.', minimum, maximum, step; FOR d IN minimum..maximum BY step LOOP INSERT INTO geodata.buffers (description, radius, geom) ( SELECT description, d, ST_Transform ( ST_Difference ( ST_Buffer(ST_Transform(geom, 3857), (d * 1000), 'quad_segs=90'), ST_Buffer(ST_Transform(geom, 3857), ((d * 1000) - (step * 1000)), 'quad_segs=90') ), 4326 )::geometry(Polygon,4326) FROM schema.tablename ); quantity := quantity + 1; END LOOP; RAISE NOTICE 'Generated % concentric ring buffers.', quantity; RETURN quantity; END; $$ LANGUAGE plpgsql;
The first line takes care of removing any older versions of this function, something that made my life a lot easier during development, since I didn’t have to manually remove it before testing a newer version.
As you can see from the header, the function takes three arguments: the minimum (starting) buffer size, the maximum (end) buffer size and the step size. All these have to be given in kilometers. You can have a look at the actual function call below to see what that looks like in real life.
The main part of the function is a for loop that steps form the minimum buffer size to the maximum buffer size in the bespoke step size. The main workhorse of this loop, and of the whole function as such, is the SQL statement that calculates the difference between two buffers: one with the current buffer size of that pass, and one with a size one step size smaller (eg. 10km – 5km). The PostGIS functions used here are
ST_Buffer to create the buffers, and
ST_Difference to calculate the difference (i.e. stamp out the hole of the donut).
In addition you can see that there are some coordinate transformations using
ST_Transform going on as well. These are necessary since
ST_Buffer does not allow for passing a unit of measurement, but will always use the distance unit of the SRS of the source data. WGS84 (EPSG 4326) uses lat/lon but I wanted to be able to create buffers in sizes specified in kilometers, so I had to temporarily transform to a metric SRS, in this case EPSG 3857). Of course after the actual buffering and differencing I had to re-transform it to WGS84 to match my other data.
Also, you might have noticed the
quad_segs=90 in the
ST_Buffer calls and wondered what that is. According to the official PostGIS documentation the
quad_segs parameter defines the “number of segments used to approximate a quarter circle” and defaults to
8. I found the resulting “circles” to be very choppy and therefore opted for a setting of
90, since that should create on circle segment per degree, which should be sufficiently precise for the kind of macro-scale analysis I wrote the function for. Feel free to change all these parameters to values that fit your needs – especially the SRS settings should match your database’s settings (see above) as do of course the schema and table names, which I decided to hardcode for simplicity.
Finally, this is how that new function can be used:
SELECT concentric_ring_buffers(5, 50, 5);
In this case I asked the function to create 10 concentric buffers with a radius difference of 5 km each between 5 km and 50 km around my point features. Below you can see the output result, put on the beautiful Stamen Toner background map (under CC-BY-3.0).
My name is Ferran Ferrer, found very useful your post.
I’m trying to adapt your code to negative buffers.
SELECT concentric_ring_buffers(0, -5, -1);
‘By value in a cycle FOr must be higher than 0’
Apart of this wants to change control of loops since resulting polygon is less than a size or similar, since I’m doing negative concentric buffers, imagine when last buffer has less area than a quantity or similar to 0, the loop stops, any suggestion?
Thanks in advance.
Modifying my request, just found some solutions like using @ step for absolute values but found that looping is more complex when substracting since operationally is:
1.Substract-st_difference- ‘buffer -1km’ to ‘original polygon’, result1
2.Substract ‘buffer -2km’ to result 1 and so on till results is null.
Using absolute values in your loop is exactly what I would have suggested in your case. But like you wrote, you’d definitely have to add a routine to the loop that checks for cases when the resulting buffer polygon would have a size of 0 or less. I’m not even sure what would be the best way to handle those cases – stop the loop, break the whole process? Probably depends on your use case. I never had this issue since I was buffering point geometries, which by design can’t have negative buffers…
[…] .Un poco de Carto para representar la Fibra y el área de influencia de las centrales en base a anillos concentricos de Konstantin Greger (futuro […]