Create buffer in meters in WGS84 layer using PyQGIS

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