I'm trying to create a model on QGIS 2.18.4 which will allow me to create multiple fixed distance buffers on a selected point(500m, 1km, 2km, 3km,...,50km) with a few clicks, rather than having to do it one by one using the MMQGIS Create Buffer function.
Because I'm working on WGS84, my solution below reprojects the selected point to World Equidistant Conic which uses meters as its distance unit. I then create the buffers, then reproject the created buffers back to WGS84.
My question is, why does my final output appear different from when I create it with MMQGIS? (See screenshot of Singapore below. Purple circle is a 5km buffer I've created using MMQGIS, and the blue circle is what I've created using my model.) I've looked at the MMQGIS python code (see code below) and it appears to apply the same workflow and CRS I've used in my model.
def mmqgis_buffer_geometry(geometry, meters):
if meters <= 0:
return None
# To approximate meaningful meter distances independent of the original CRS,
# the geometry is transformed to an azimuthal equidistant projection
# with the center of the polygon as the origin. After buffer creation,
# the buffer is transformed to WGS 84 and returned. While this may introduce
# some deviation from the original CRS, buffering is assumed in practice
# to be a fairly inexact operation that can tolerate such deviation
wgs84 = QgsCoordinateReferenceSystem()
wgs84.createFromProj4("+proj=longlat +datum=WGS84 +no_defs")
latitude = str(geometry.centroid().asPoint().y())
longitude = str(geometry.centroid().asPoint().x())
#proj4 = "+proj=aeqd +lat_0=" + str(geometry.centroid().asPoint().y()) + \
# " +lon_0=" + str(geometry.centroid().asPoint().x()) + \
# " +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
# For some reason, Azimuthal Equidistant transformation noticed to not be
# working on 10 July 2014. World Equidistant Conic works, but there may be errors.
proj4 = "+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=60 +lat_2=60 " + \
"+x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
azimuthal_equidistant = QgsCoordinateReferenceSystem()
azimuthal_equidistant.createFromProj4(proj4)
transform = QgsCoordinateTransform(wgs84, azimuthal_equidistant)
geometry.transform(transform)
newgeometry = geometry.buffer(meters, 7)
wgs84 = QgsCoordinateReferenceSystem()
wgs84.createFromProj4("+proj=longlat +datum=WGS84 +no_defs")
transform = QgsCoordinateTransform(azimuthal_equidistant, wgs84)
newgeometry.transform(transform)
return newgeometry
def mmqgis_buffer_point(point, meters, edges, rotation_degrees):
if (meters <= 0) or (edges < 3):
return None
# Points are treated separately from other geometries so that discrete
# edges can be supplied for non-circular buffers that are not supported
# by the QgsGeometry.buffer() function
wgs84 = QgsCoordinateReferenceSystem()
wgs84.createFromProj4("+proj=longlat +datum=WGS84 +no_defs")
# print "Point " + unicode(point.x()) + ", " + unicode(point.y()) + " meters " + unicode(meters)
polyline = []
for edge in range(0, edges + 1):
degrees = ((float(edge) * 360.0 / float(edges)) + rotation_degrees) % 360
polyline.append(mmqgis_endpoint(point, meters, degrees))
return QgsGeometry.fromPolygon([polyline])
Why is this happening? What can I do to get the same result as when using MMQGIS?