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.

Improve this question

1 Answer 1

Reset to default
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})
Improve this answer
5
  • 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
    Commented Feb 14, 2020 at 8:38
  • I edited my answer, I give you another code, try it and give me a feedback.
    – Vincent Bré
    Commented Feb 14, 2020 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é
    Commented Feb 14, 2020 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
    Commented Feb 14, 2020 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é
    Commented Feb 14, 2020 at 10:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.