When you want to buffer features that are spread across a large area (such as global layers), there is no suitable projection that can give you accurate results. This is the classic case for needing Geodesic Buffers - where the distances are measured on an ellipsoid or spherical globe. This post explains the basics of geodesic vs. planar buffers well.

QGIS lacks a way to do geodesic buffers natively. But one can approximate them by using a custom Azimuthal Equidistant projection for each point. Azimuthal Equidistant projection has the useful property that all points on the map are at proportionally correct distances from the center point. So we can write a custom processing script that implements the following algorithm for very accurate buffers - even for global layers.

  • For each feature in a layer, create a custom azimuthal equidistant projection centered at that geometry.
  • Transform the geometry to the custom projection and do a planar buffer.
  • Reverse transform the buffer to the original projection
  • Repeat for each feature

MMQGIS plugin has a similar implementation in the buffer tool, but the current version uses the World Equidistant Conic projection (instead of custom azimuthal equidistant projection for each feature) which is not as accurate.

This post explains how to write a QGIS Processing script that implements this logic. It also serves as a reference for QgsProcessingFeatureBasedAlgorithmclass which is the preferred and simplified way to implement algorithm that work on each feature independently. If you just want to use the script and are not interested in the implementation details, you can do the following and skip rest of the post

  • Download the equidistance_buffer.py script
  • In QGIS, open Settings → User Profiles → Open Active Profile Folder
  • Copy the equidistance_buffers.py script to processing → scripts folder.
  • Restart QGIS and launch the script from Processing Toolbox → Scripts → scripts → Equidistance Buffer

Tip: While buffering global layers, sometimes you may find that your buffer polygon is crossing the dateline, causing artifacts when displaying them in some projections. You can avoid this by splitting the buffer polygon into 2 halves at the dateline. ogr2ogr tool can do this automativcally using the -wrapdateline option

Processing Script

Writing Processing Python scirpts in QGIS involves using the QgsProcessingAlgorithm class (see my tutorial if you are looking for an introduction). For algorithms which operate on features one-by-one, there’s a class *QgsProcessingFeatureBasedAlgorithm,*which makes for a much simpler implementation. Since our problem operates on each feature independently, we can use this class instead.

QgsProcessingFeatureBasedAlgorithmrequires a few different methods to be implemented for initializing, preparing and processing of the algorithm. We will look into each of them

Inputs and Outputs. There is no need to declare the input source or output sink, as these are automatically created by the class. You just need to declare what would be your output layer type, so an appropriate output sink can be created. Similarly, you can define valid input source types to filter the selection of layers that the user can choose from. Here our input can any vector layer and the output is a polygon layer.

def inputLayerTypes(self):
         return [QgsProcessing.TypeVector]
def outputWkbType(self, input_wkb_type):
    return QgsWkbTypes.Polygon

If you need any extra parameters apart from input and output, you need to define them in the initParametersmethod. For our script, we need to take a distance input from the user.

 def initParameters(self, config=None):
        self.addParameter(
            QgsProcessingParameterDistance(
                self.DISTANCE,
                'Buffer Distance (meters)',
                defaultValue=10000
                ))

We need to store the user supplied distance value in an instance variable (referred by self.distance) so it can be accessed when processing each feature. This is done in the prepareAlgorithmmethod which is for evaluating parameter values. We can also check if the chosen input layer is in a Geographic CRS and raise an error if it is not.

def prepareAlgorithm(self, parameters, context, feedback):
        self.distance = self.parameterAsDouble(
            parameters,
            self.DISTANCE,
            context)
        source = self.parameterAsSource(parameters, 'INPUT', context)
        self.source_crs = source.sourceCrs()
        if not self.source_crs.isGeographic():
            feedback.reportError('Layer CRS must be a Geograhpic CRS for this algorithm')
            return False
        return super().prepareAlgorithm(parameters, context, feedback)

Next is the implementation of the processFeaturemethod. This method should contain the main logic of the script. You just need to define how a single feature should be processed, and QGIS takes care of iterating and writing the results to the output. Note that the processFeaturemethod can create multiple output features, so the return value is a list. Here we implement the logic of transforming the feature geomety to a azimuthal equidistant projection, buffering it and converting the resulting geometry to the original projection.

  def processFeature(self, feature, context, feedback):
        geometry = feature.geometry()
        # For point features, centroid() returns the point itself
        centroid = geometry.centroid()
        x = centroid.asPoint().x()
        y = centroid.asPoint().y()
        proj_string = 'PROJ4:+proj=aeqd +ellps=WGS84 +lat_0={} +lon_0={} +x_0=0 +y_0=0'.format(y, x)
        dest_crs = QgsCoordinateReferenceSystem(proj_string)
        xform = QgsCoordinateTransform(self.source_crs, dest_crs, QgsProject.instance())
        geometry.transform(xform)
        buffer = geometry.buffer(self.distance, self.SEGMENTS, self.END_CAP_STYLE, self.JOIN_STYLE, self.MITER_LIMIT)
        buffer.transform(xform, QgsCoordinateTransform.ReverseTransform)
        feature.setGeometry(buffer)
        return [feature]

Below is the complete script for reference.

from PyQt5.QtCore import QCoreApplication
from qgis.core import (QgsProcessing, QgsProcessingFeatureBasedAlgorithm,
    QgsProcessingParameterDistance, QgsPoint, QgsFeature, QgsGeometry,
    QgsWkbTypes, QgsCoordinateReferenceSystem, QgsCoordinateTransform, QgsProject)

class EquidistanceBuffer(QgsProcessingFeatureBasedAlgorithm ):
    """
    This algorithm takes a vectorlayer and creates equidistance bufers.
    """
    DISTANCE = 'DISTANCE'
    SEGMENTS = 5
    END_CAP_STYLE = QgsGeometry.CapRound
    JOIN_STYLE = QgsGeometry.JoinStyleRound
    MITER_LIMIT = 2.0

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)
    def createInstance(self):
        return EquidistanceBuffer()
    def name(self):
        return 'equidistance_buffer'
    def displayName(self):
        return self.tr('Equidistance Buffer')
    def group(self):
        return self.tr('scripts')
    def groupId(self):
        return 'scripts'
    def outputName(self):
        return self.tr('buffers')

    def shortHelpString(self):
        return self.tr(
            'This algorithm creates equidistant buffers for vector layers. \
            Each geometry is transformed to a Azimuthal Equidistant projection \
            centered at that geometry, buffered and transformed back to the \
            original projection.')
    def __init__(self):
        super().__init__()

    def inputLayerTypes(self):
        return [QgsProcessing.TypeVector]

    def outputWkbType(self, input_wkb_type):
        return QgsWkbTypes.Polygon

    def initParameters(self, config=None):
        self.addParameter(
            QgsProcessingParameterDistance(
                self.DISTANCE,
                'Buffer Distance (meters)',
                defaultValue=10000
                ))

    def prepareAlgorithm(self, parameters, context, feedback):
        self.distance = self.parameterAsDouble(
            parameters,
            self.DISTANCE,
            context)
        source = self.parameterAsSource(parameters, 'INPUT', context)
        self.source_crs = source.sourceCrs()
        if not self.source_crs.isGeographic():
            feedback.reportError('Layer CRS must be a Geograhpic CRS for this algorithm')
            return False
        return super().prepareAlgorithm(parameters, context, feedback)
    def processFeature(self, feature, context, feedback):
        geometry = feature.geometry()
        # For point features, centroid() returns the point itself
        centroid = geometry.centroid()
        x = centroid.asPoint().x()
        y = centroid.asPoint().y()
        proj_string = 'PROJ4:+proj=aeqd +ellps=WGS84 +lat_0={} +lon_0={} +x_0=0 +y_0=0'.format(y, x)
        dest_crs = QgsCoordinateReferenceSystem(proj_string)
        xform = QgsCoordinateTransform(self.source_crs, dest_crs, QgsProject.instance())
        geometry.transform(xform)
        buffer = geometry.buffer(self.distance, self.SEGMENTS, self.END_CAP_STYLE, self.JOIN_STYLE, self.MITER_LIMIT)
        buffer.transform(xform, QgsCoordinateTransform.ReverseTransform)
        feature.setGeometry(buffer)
        return [feature]