5

In my QGIS 3 plugin I have this code that returns the list of all features intersected with a buffer calculated from two input values: the ID of the feature in the center of the buffer and a radius length. The problem is that the layer is in CRS WGS84, that uses degrees as unit. I need to calculate the buffer using meters as unit.

Is there a way to calculate this buffer in meters without change layer's CRS?

expr = QgsExpression("id = '"+InputID+"'")
it = layer.getFeatures(QgsFeatureRequest(expr))
ids = [i.id() for i in it]
layer.selectByIds(ids)
selection = layer.selectedFeatures()
for feat in selection:
    geom = feat.geometry().centroid()
    buff = geom.buffer(inputDistance,5)
intersectedFeaturesList=[]
for feature in layer.getFeatures():
    if feature.geometry().intersects(buff):
        intersectedFeaturesList.append(str(feature.attributes()[0])) 
return intersectedFeaturesList

EDIT: I solved in this way, but I'm open to new suggestions.

import processing
from processing.core.Processing import Processing

Processing.initialize()
parameter = {'INPUT': layer, 'TARGET_CRS': 'EPSG:3857','OUTPUT': 'memory:Reprojected'}
result = processing.run('native:reprojectlayer', parameter)
layer3857 = result['OUTPUT']
expr = QgsExpression("id = '"+InputID+"'")
it = layer3857.getFeatures(QgsFeatureRequest(expr))
ids = [i.id() for i in it]
layer3857.selectByIds(ids)
selection = layer3857.selectedFeatures()
for feat in selection:
    geom = feat.geometry().centroid()
    buff = geom.buffer(inputDistance,5)
intersectedFeaturesList=[]
for feature in layer3857.getFeatures():
    if feature.geometry().intersects(buff):
        intersectedFeaturesList.append(str(feature.attributes()[0])) 
return intersectedFeaturesList
1

Here is an alternative approach which avoids using processing- use the QgsCoordinateTransform class to apply a CRS transformation to feature geometries prior to creating buffers and testing for intersections.

Disclaimer* Although its units are in meters, just be aware that the Web Mercator (EPSG:3857) CRS is notoriously inaccurate for distance/ area calculations, more so further from the equator - so use it at your discretion.

# imports for use in a plugin
from qgis.core import QgsCoordinateReferenceSystem, QgsCoordinateTransform

# variables created for testing
input_id = 1
input_distance = 50000 # 50km
layer = QgsProject().instance().mapLayersByName('Your_Layer_Name')[0]

layer.selectByIds([input_id])
all_feats = [f for f in layer.getFeatures()]

# in case you want to exclude input feature from intersection test against...
#...its own buffered centroid use line below
#all_feats = [f for f in layer.getFeatures() if f.id() != input_id]

input_feat = layer.getFeature(input_id)
# construct a QgsCoordinateReferenceSystem object
new_crs = QgsCoordinateReferenceSystem('EPSG:3857')
# construct a QgsCoordinateTransform object
xtransform = QgsCoordinateTransform(layer.crs(), new_crs, QgsProject().instance())
geom = input_feat.geometry()
# transform input feature geometry
geom.transform(xtransform)

intersecting_features = []

for f in all_feats:
    f_geom = f.geometry()
    # transform geometry of each feature in current iteration before...
    # ...testing for intersection
    f_geom.transform(xtransform)
    if f_geom.intersects(geom.centroid().buffer(input_distance, 5)):
        intersecting_features.append(f)

print(intersecting_features)
3
  • Hi thanks again for your alternative solution and for the suggestion about EPSG:3857. Which CRS would you use?
    – Lorenzo
    Feb 8 '20 at 11:22
  • 3
    EPSG:3857 is a hack, and may generate distances which are woefully inaccurate. Assuming your buffer is a couple of metres or kilometres, the differences in output from "real" projections will be minimal. The right answer is most likely whatever is used by cadastral land surveyors in your country if the application is for property/land issues, and whatever is used for official topo maps otherwise. If not worth the complexity deciding, I'd use the UTM zone for your country picked from here: downloads.gvsig.org/download/geodata/EPSG_codes/…
    – Houska
    Feb 8 '20 at 11:42
  • 1
    Agree with what @Houska said. Selecting a CRS is not a trivial topic and it depends on factors including your region and the size of your area of interest etc. I don't purport to have the greatest expertise in this area, but I think that a projected UTM CRS for the appropriate zone as suggested by Houska would probably be a suitable choice.
    – Ben W
    Feb 8 '20 at 12:11
3

An option is to make use of the geographiclib package to define a function that creates a "geodesic buffer".

This function solves the Direct problem of geodesy to find the coordinates of each point at a defined ellipsoidal distance from a coordinates pair, for initial azimuthes defined by the number of segments that the buffer polygon geometry must have (for each quart circle).

The function is not handling errors, so please send valid arguments to it.

From the console, a working code may looks like the following:

# requires geographiclib 
# (https://geographiclib.sourceforge.io/html/python/index.html)
from geographiclib.geodesic import Geodesic

def geodesic_buffer(lon,lat,distance,segments):
    """Returns a Polygon Geometry,
    geodesic (WGS84) buffer from a pair of coordinates.

    lon --- longitude (WGS84) of a point
    lat --- latitude (WGS84) of a point
    distance --- ellipsoidal distance for the buffer
    segments --- segments (each 90 degrees) of the polygon
    """
    geod = Geodesic.WGS84 # geodesic object
    angle = 90 / segments # internal angles of the polygon
    coords = []
    for i in range(segments * 4 + 1):
        direct = geod.Direct(lat,lon,i*angle,distance) #solve the Direct problem
        coords.append((direct['lon2'],direct['lat2']))
    # create the polygon geometry from the coords list
    geo_buffer = QgsGeometry.fromPolygonXY([[QgsPointXY(pair[0],pair[1]) for pair in coords]])
    # print(geo_buffer.asWkt())
    return geo_buffer


layer_name = 'test'
input_id = 1
input_distance = 50000 # meters
layer = QgsProject().instance().mapLayersByName(layer_name)[0]
feat = layer.getFeature(input_id)
geom = feat.geometry().centroid()
buff = geodesic_buffer(geom.asPoint().x(),geom.asPoint().y(),input_distance,5)
intersected_features_list=[]
for feature in layer.getFeatures():
    if feature.geometry().intersects(buff):
        intersected_features_list.append(str(feature.attributes()[0]))

print(intersected_features_list)

Another option is to solve the Direct problem with pyproj:

from pyproj import Geod

def geodesic_buffer(lon,lat,distance,segments):
    """Returns a Polygon Geometry,
    geodesic (WGS84) buffer from a pair of coordinates.

    lon --- longitude (WGS84) of a point
    lat --- latitude (WGS84) of a point
    distance --- ellipsoidal distance for the buffer
    segments --- segments (each 90 degrees) of the polygon
    """
    g = Geod(ellps='WGS84') # geodesic object
    angle = 90 / segments # internal angles of the polygon
    coords = []
    for i in range(segments * 4 + 1):
        direct = g.fwd(lon,lat,i*angle,distance) #solve the Direct problem
        coords.append((direct[0],direct[1]))
    # create the polygon geometry from the coords list
    geo_buffer = QgsGeometry.fromPolygonXY([[QgsPointXY(pair[0],pair[1]) for pair in coords]])
    # print(geo_buffer.asWkt())
    return geo_buffer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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