CODEWITHSUNDEEP
pythonpolygoncomputational-geometryshapefileqgis
Shubham Goyal
Shubham Goyal

Reputation: 1322

QGis: Finding the polygon inside which a point lies

I have a layer with polygonal features. Each feature has attributes and values. I also have a list of co-ordinates and I would like to know which feature (or polygon) the co-ordinates lie in.

Could someone please guide me on how to go about this? Is there a function in the API that can help me achieve my goal or should I use some computational geometry algorithm to do it myself? I know how to do the latter but it would save me some time if there was a built in function already.

Thanks.

Upvotes: 2

Views: 4305

Answers (2)

user2716238
user2716238

Reputation: 21

while provider.nextFeature(feature):
    if (feature.geometry().contains(QgsGeometry.fromPoint(QgsPoint(lon, lat)))):
        print 'Contained in feature %d' % feature.id()

Upvotes: 2

Shubham Goyal
Shubham Goyal

Reputation: 1322

I eventually managed to do it myself.

import sys
import os
from qgis.core import *
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches

LATITUDE = 1.29306
LONGITUDE = 103.856

QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()

layer=QgsVectorLayer("/home/shubham/SMART/zones/mtz1092p.shp", "mtz1092p", "ogr")
if not layer.isValid():
    print "Layer failed to load!"
provider = layer.dataProvider()

def findFeatureId(point):
    feat = QgsFeature()
    allAttrs = provider.attributeIndexes()
    provider.select(allAttrs)
    while provider.nextFeature(feat):
        geom = feat.geometry()
        x = geom.asPolygon()
        if len(x) == 0:
            print "Feature ID %d has no ring" % feat.id()
        else:
            codes = []
            codes.append(Path.MOVETO)
            for i in range (0, len(x[0]) - 2):
                codes.append(Path.LINETO)
            codes.append(Path.CLOSEPOLY)
            path = Path(x[0], codes)
            if (path.contains_point(point, None, 0.0)):
                print "Point contained in feature ID %d" %feat.id()

if __name__ == "__main__":
    crsSrc = QgsCoordinateReferenceSystem(4326) # WGS84
    crsDest = QgsCoordinateReferenceSystem(3414)# SVY21
    xform = QgsCoordinateTransform(crsSrc, crsDest)
    pt = xform.transform(QgsPoint(LONGITUDE, LATITUDE))
    findFeatureId(pt)

Upvotes: 1

Related Questions