Saturday 19 October 2019

pyqgis - How to create rectangular buffers around points in QGIS with Python?


Using the following code from a similar question (How to create square buffers around points in QGIS with Python?) I am able to create square buffers (5x5 degree cells) around my points.


However, I also need to make 5x10 degree cells and 10x20, is there a way I can modify the following code to get this result?


Code:


layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]


epsg = layer.crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=x:real&field=y:real&field=point_id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
'square_buffer',
'memory')

prov = mem_layer.dataProvider()


for i, feat in enumerate(feats):
point = feat.geometry().asPoint()
new_feat = QgsFeature()
new_feat.setAttributes([i, point[0], point[1], feat.id()])
tmp_feat = feat.geometry().buffer(1000, -1).boundingBox().asWktPolygon()
new_feat.setGeometry(QgsGeometry.fromWkt(tmp_feat))
prov.addFeatures([new_feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)


Answer



If the code is modified as follow, it can be obtained rectangular buffers as you expected.


layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

epsg = layer.crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=x:real&field=y:real&field=point_id:integer""&index=yes"


mem_layer = QgsVectorLayer(uri,
'rectangular_buffer',
'memory')

prov = mem_layer.dataProvider()

for i, feat in enumerate(feats):
point = feat.geometry().asPoint()
new_feat = QgsFeature()
new_feat.setAttributes([i, point[0], point[1], feat.id()])

bbox = feat.geometry().buffer(1000, -1).boundingBox()
tmp_feat = bbox.asWktPolygon()
xmin1,ymin1,xmax1,ymax1 = bbox.toRectF().getCoords()
xmin2,ymin2,xmax2,ymax2 = feat.geometry().buffer(2000, -1).boundingBox().toRectF().getCoords()
p1 = QgsPoint(xmin1, ymax2)
p2 = QgsPoint(xmax1, ymin2)
new_ext = QgsRectangle(p1,p2)
new_tmp_feat = new_ext.asWktPolygon()
new_feat.setGeometry(QgsGeometry.fromWkt(new_tmp_feat))
prov.addFeatures([new_feat])


QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

After running the code at the Python Console of QGIS, I got:


enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...