How to create variable distance buffers around polygons in Python

Question

I have a shapefile of polygons and I would like to create buffers around the edges of each polygon. However, rather than use a single buffer distance for all of the polygons, I am trying to find a way to create variable distance buffers. Specifically, within the attribute table for this layer, I have a column of Distances where each polygon has assigned to it its own unique distance. I want to create buffers around each polygon that use each polygon's unique distance value as its unique buffer distance.

I have found so far from "Python GDAL/OGR Cookbook 1.0 documentation" (linked here: http://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html?highlight=buffer) is that I can use the following to create buffers:

createBuffer(inputfn, outputBufferfn, bufferDist):

However, this entails using the same buffer distance for each polygon. What I am trying to figure out specifically is, how can I use this createBuffer() function (or another function if this one is not appropriate here) to use my Distances column to provide a unique buffer distance for each polygon?

• jbalk: On the cookbook page you linked, it defines the buffer function. I would modify that function for your case. Remove bufferDist parameter. Then, in the for feature in inputLyr:... loop, add something like bufferDist = feature.getField("Distances").

Answer 1

Modified version of createBuffer function using a column as the distance value is here:

import ogr, os
def createBuffer(inputfn, outputBufferfn, bufferField):
 inputds = ogr.Open(inputfn)
 inputlyr = inputds.GetLayer()
 shpdriver = ogr.GetDriverByName('ESRI Shapefile')
 if os.path.exists(outputBufferfn):
 shpdriver.DeleteDataSource(outputBufferfn)
 outputBufferds = shpdriver.CreateDataSource(outputBufferfn)
 bufferlyr = outputBufferds.CreateLayer(outputBufferfn, geom_type=ogr.wkbPolygon)
 featureDefn = bufferlyr.GetLayerDefn()
 for feature in inputlyr:
 ingeom = feature.GetGeometryRef()
 bufferValue = feature.GetField(bufferField)
 geomBuffer = ingeom.Buffer(bufferValue)
 outFeature = ogr.Feature(featureDefn)
 outFeature.SetGeometry(geomBuffer)
 bufferlyr.CreateFeature(outFeature)
 outFeature = None
inputfn = 'input.shp'
outputBufferfn = 'output.shp'
bufferField = 'Distances'
createBuffer(inputfn, outputBufferfn, bufferField)

• LostinSpatialAnalysis: I am testing this out now, and the function itself is created fine. However, when I enter inputfn = 'input.shp'', 'outputBufferfn = 'out.shp', 'bufferField = 'Distances'', 'createBuffer(inputfn, outputBufferfn, bufferField)' I get the following error: RuntimeError: not a string. I checked an Distances` is listed as float64 via inputfn.dtypes. Could this mean that there is still something wrong with the column I am referencing for bufferField? Also, rather than using input.shp for inputfn I am just referencing a geopandas dataframe.
• LostinSpatialAnalysis: When I do use 'inputfn = 'input.shp'' I get the error: AttributeError: 'NoneType' object has no attribute 'CreateLayer'

Answer 2

You can use the GeoPandas library to read in your Shapefile as a GeoDataFrame.

Then you can make use of the GeoDataFrame's vectorized buffer function, as suggested in this answer by Encomium.

Here's a quick mock-up of what it would look like.

import geopandas as gpd
# Reading in the input shapefile.
input_df = gpd.read_file('path/to/my_shapefile.shp')
# Reading the input CRS. 
input_crs = input_df.crs
# Reprojecting the data. If needed, substitute "WXYZ" with 
# relevant EPSG code. Learn more here: https://spatialreference.org/ref/epsg/
# Comment out if unnecessary for your specific case.
input_df = input_df.to_crs('epsg:WXYZ')
# Creating the variable-sized buffer
input_df['buffer'] = input_df.buffer(input_df['Distances'])
# Dropping the original geometry and setting the new geometry
buff_df = input_df.drop(columns=['geometry']).set_geometry('buffer')
# Reprojecting the buffered data back to the CRS used in 
# the original input shapefile
buff_df = buff_df.to_crs(input_crs)
# Exporting the output to a shapefile
buff_df.to_file('path/to/buffered_output.shp')

This will yield a new shapefile with the same number of rows/features, but each one will be buffered using the value from the Distances column.

Important note: Mind the CRS!

GeoPandas does not know how to measure distances across different projection/coordinate systems. Therefore, you should pay close attention to the compatibility between the CRS system being used in your input shapefile and the distance units being used in your Distances column. Stated clearly, you need to make sure that both the CRS and the Distances column have the same units.

For example, if the input shapefile has coordinates coded in WGS84/EPSG4326 (i.e., decimal degree coordinates on an ellipsoid), and your Distances column has values in some real-world length unit (for example, feet), the results will not look like you expect. You first will need to convert the data from the input shapefile into a projected CRS.

Thankfully, this is quite easy in GeoPandas. You can just use input_df = input_df.to_crs('epsg:WXYZ'), where you replace the WXYZ with the EPSG code of a projected CRS that works for your specific case. You can learn more about EPSG codes here and here.

Edit

My original answer used a less efficient method. The current answer uses a vectorized approach suggested here by @Encomium.

• Paul H: Links might go away. I think you should include (with attribution) the definition of the make_variable_buffer function in this answer
• Felipe D.: Good point! Thanks! I've added the definition to the body of the main answer. Also, the attribution is easy: I'm the one who wrote the other answer too! Hahahaha...
• Felipe D.: Actually, scratch that. Looking at the original question I linked to, I found an even better approach. I've edited my answer and given the proper attribution to the original solution.

Answer 3

This is something where SQL queries work nicely. Haven't tested it but this might be close to what you want: select *, st_Buffer(geometry, my_vect.Distance) from my_vect

You can do this using ogr2ogr on the command line:

ogr2ogr \
 -sql "select *, ST_Buffer(geometry, my_vect.Distance) from my_vect" `\
 -dialect SQLite output_vect.shp my_vect.shp

If you want to use it from within Python, gdal.VectorTranslate will allow you to translate (sorry! ;-D) the previous shell command to a Python command.

Note that Distance needs to be in the same units as the CRS of the geometry you are expanding, but you can read more about this on the documentation