8

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.

4 Answers 4

6

Few features

What you probably want to do is:

  1. Create a list of QgsPoint from your coordinates
  2. Iterate over all your layer features (polygons)
  3. Loop over the list of points (within the iteration)
  4. Call feature.geometry().contains( point ) to check if you've got a match

Of course you can now improve performance if you e.g. know, that a point can only be contained by one polygon, you can remove a point from the list, once the appropriate polygon has been found.

Lots of features (Using SpatialIndex)

As pointed out in the comments, a spatial index can be used to speed up the process significantly.

The steps here would be

  1. Create a list of QgsPoint from your coordinates
  2. Create a QgsSpatialIndex
  3. Iterate over all your layer features (polygons)
  4. Add the features to your index with insertFeature
  5. Iterate over all your points
  6. Call index.intersects( QgsRectangle( point, point ) ) to check if you've got a match

There is also a code example by NathanW

4
  • Ah, I didn't know about feature.geometry.contains(point) call. I used mathplotlib. pastebin.com/61LSeMWv Please do forgive the untidiness of my code. I am in a hurry and will clean it later. May 30, 2013 at 8:47
  • I haven't implemented your solution and so can't say for certain whether it works or not. But I do believe it should and so I am marking this as the correct answer :) May 30, 2013 at 8:51
  • 2
    Could this be sped up using a QgsSpatialIndex? May 30, 2013 at 9:06
  • 1
    I would highly recommend using a QgsSpatailIndex. There is an example at nathanw.net/2013/01/04/…
    – Nathan W
    May 31, 2013 at 5:41
6

First of all you have to import the list of coordinates into your project. This tutorial explains well how to do that: http://qgis.spatialthoughts.com/2012/01/importing-spreadsheets-or-csv-files-to.html

When you have both layers (polygones and points) into your project, go to vector > data management tools > join attributes by location

enter image description here

You get a window where you can define which layers you want to combine:

enter image description here

  • Set your point layer als the 'target vector layer'.
  • Set your polygon layer as 'join vector layer'.
  • Define your Output Shapefile (there will be a new one created. So if you missed, the original data is preserved).
  • You can choose to keep all the data in the new shapefile, even if there is no match with a polygon: check 'keep all records (including non-matching target records)'

Click 'OK'. The new shapefile is created and you will be asked 'Would you like to add the new layer to the TOC?' Click again OK.

Open the attributetable of the new added shapefile and you will see that all the features of the corresponding polygon are added to the point that lays into that polygon.

2

A simpler way to do this using PyQGIS. I figured that you can construct a QgsRectangle object with a single point and use it with QgsFeatureRequest to filter features from the layer that intersect it.

point = QgsPoint(10, 10)
# Construct a QgsRectange with the same point
rect = QgsRectangle(point, point)
req = QgsFeatureRequest()
req.setFilterRect(rect)
# You get the feature that intersects the point
f = layer.getFeatures(req).next()
0

In QuantumGIS you could add the list of coordinates with the function 'add delimited text layer' (if it is a csv file). Also add the polygons. Then you can do an 'intersect' or 'points in polygon'.

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.