Selecting features inside a polygon - PyQGIS

Question

I have just started working with PyQGIS and I would like to ask this question. I have a rather big dataset (cell phone towers) and I want to narrow it in the towers that exist inside a certain country and I want to do it with PyQGIS. I used that piece of code:

lyrPts = QgsVectorLayer("/path/Cell-tower.shp", "CellTowers", "ogr")

lyrPoly =QgsVectorLayer("/path/GRC_adm0.shp", "GRC", "ogr")

QgsProject.instance().addMapLayers([lyrPts,lyrPoly])

ftsPoly = lyrPoly.getFeatures()

for feat in ftsPoly:
    geomPoly = feat.geometry()
    bbox = geomPoly.boundingBox()
    req = QgsFeatureRequest()
    filterRect = req.setFilterRect(bbox)
    featsPnt = lyrPts.getFeatures(filterRect)
    for featPnt in featsPnt:
        if featPnt.geometry().within(geomPoly):
            lyrPts.select(featPnt.id())

iface.setActiveLayer(lyrPoly)
iface.zoomToActiveLayer()

Although when I run it I don't get any errors I also don't get any results even when I left it running all night.

Can anyone help me and tell me if the code is good?

For reference the cell tower data I got them from mylnikov.org as a .csv file and made it to .shp file and the polygon from diva-gis.org.

Vincent Bré · Answer 1 · 2020-02-14T09:47:43.610

2

I replaced the first two lines of your code because my layers are already in the QGIS project and I removed the third one. I execute the following code :

lyrPts = QgsProject.instance().mapLayersByName('point')[0]
lyrPoly = QgsProject.instance().mapLayersByName('poly')[0]

ftsPoly = lyrPoly.getFeatures()

for feat in ftsPoly:
    geomPoly = feat.geometry()
    bbox = geomPoly.boundingBox()
    req = QgsFeatureRequest()
    filterRect = req.setFilterRect(bbox)
    featsPnt = lyrPts.getFeatures(filterRect)
    for featPnt in featsPnt:
        if featPnt.geometry().within(geomPoly):
            lyrPts.select(featPnt.id())

iface.setActiveLayer(lyrPoly)
iface.zoomToActiveLayer()

I get the result instantly or in few seconds. The result is the following :

I give you another code thats run faster which use the tool of QGIS GUI:

lyrPts = QgsProject.instance().mapLayersByName('point')[0] # Point layer
lyrPoly = QgsProject.instance().mapLayersByName('poly')[0] # Polygon layer
processing.run("native:selectbylocation", {'INPUT':lyrPts,'PREDICATE':[6],'INTERSECT':lyrPoly,'METHOD':0})

edited Feb 14 '20 at 09:47

answered Feb 14 '20 at 08:08

Vincent Bré

4,090
7
23

Yes they are in the same projection. The thing is that the cell tower dataset is huge. What I am trying to do is to extract the cell towers by country. So I need all the cell-towers inside the polygon representing the country. When I run the code I don't get an error but it takes way to much time. When I do it in Qgis with the GUI it doesn't take so much time. So I was wandering if there was a problem with the code making it go so slow.
Anyway. Thank you very much for your help
– ovidius Feb 14 '20 at 08:38
I edited my answer, I give you another code, try it and give me a feedback. – Vincent Bré Feb 14 '20 at 08:53
You get an error ? Are you sure some of your points are within your polygons? You're talking about one big data set, how many features? SQL may seem like a good alternative. – Vincent Bré Feb 14 '20 at 09:50
It works. Thank you very much. Would it be to impolite to ask some help in how to present it on a new layer or and export it to a .shp file again using code? – ovidius Feb 14 '20 at 10:01
You're welcome, you can accept my answer. For your question, you can follow this link : https://gis.stackexchange.com/questions/80292/how-can-i-create-a-vector-layer-from-selected-features-with-pyqgis/262405. Good luck! – Vincent Bré Feb 14 '20 at 10:04

Selecting features inside a polygon - PyQGIS

1 Answers1