2

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.

New contributor
ovidius is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
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 :

enter image description here

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})
  • 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 at 8:38
  • I edited my answer, I give you another code, try it and give me a feedback. – Vincent Bré Feb 14 at 8: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 at 9: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 at 10:01
  • You're welcome, you can accept my answer. For your question, you can follow this link : gis.stackexchange.com/questions/80292/…. Good luck! – Vincent Bré Feb 14 at 10:04

Your Answer

ovidius is a new contributor. Be nice, and check out our Code of Conduct.

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.