1

Using a spatial query in QGIS Desktop to select features of layer MTA_ACQ_TE_polygon that are within features of layer nts_snrc_50k works, but I am attempting to do this programatically in PyQGIS, with no luck.

Screen shot of QGIS Desktop

And this is my attempt:

from qgis.core import *

app=QgsApplication([], False)
app.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
app.initQgis()

mapsheet = QgsVectorLayer("E:/MTOdata/archives/0/MTA_ACQ_TE/nts_snrc_50k.shp", \
    "mapsheet", "ogr")
tenures = QgsVectorLayer("E:/MTOdata/archives/0/MTA_ACQ_TE/MTA_ACQ_TE_polygon.shp", \
    "tenures", "ogr")

if not mapsheet.isValid() or not tenures.isValid():
    print("failed to open shapefiles")
    exit(0)

print("-> loaded mapsheet (%s), tenures (%s)" % (mapsheet.crs().authid(), tenures.crs().authid()))
print("-> finding mapsheet")

for m in mapsheet.getFeatures():
    if m["NTS_SNRC"] == "094E06":
        print("\tfound it!")
        sheet = m.geometry()
        break

print("-> finding claims within")

for t in tenures.getFeatures():
    tenure = t.geometry()
    if tenure.within(sheet):
        print("\tgot one: %s" % t["OBJECTID"])

print("-> done")

Which provides the output:

-> loaded mapsheet (EPSG:4140), tenures (EPSG:3005)
-> finding mapsheet
        found it!
-> finding claims within
-> done

It fails to find any features of MTA_ACQ_TE_polygon within nts_snrc_50k. Why?

  • Have you tried contains() instead of within() (you'll have to switch tenure and sheet)? – ArMoraer Apr 7 '16 at 19:54
  • Yes I have tried variants of contains, within and even intersect in an effort to debug. – Jack Morton Apr 7 '16 at 19:56
  • And have you tried to reproject one of the layers into the same CRS as the other? – ArMoraer Apr 7 '16 at 20:03
1

I followed ArMoraer's advice, and have reprojected each instance of tenure geometry to match the projection of nts_snrc_50k.

The code looks like this:

alb_crs = QgsCoordinateReferenceSystem(3005, QgsCoordinateReferenceSystem.EpsgCrsId)
longlat_crs = QgsCoordinateReferenceSystem(4140, QgsCoordinateReferenceSystem.EpsgCrsId)
tr = QgsCoordinateTransform(alb_crs, longlat_crs)

And like this:

for t in tenures.getFeatures():
    t.geometry().transform(tr)
    tenure = t.geometry()
    if tenure.within(sheet):
        print("\tgot one: %s" % t["OBJECTID"])

Your Answer

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.