官术网_书友最值得收藏!

Filtering a layer by geometry

In this recipe, we'll perform a spatial operation to select a subset of a point layer based on the points contained in an overlapping polygon layer. We'll use shapefiles in both cases, with one being a point layer and the other a polygon. This kind of subset is one of the most common GIS operations.

Getting ready

We will need two new shapefiles that have not been used in previous recipes. You can download the point layer from https://geospatialpython.googlecode.com/files/MSCities_Geo_Pts.zip.

Similarly, you can download the geometry layer from https://geospatialpython.googlecode.com/files/GIS_CensusTract.zip.

Unzip these shapefiles and place them in a directory named ms within your qgis_data directory, within your root or home directory.

How to do it...

In this recipe, we will perform several steps to select features in the point layer that fall within the polygon layer, as follows:

  1. First, load the point layer:
    lyrPts = QgsVectorLayer("/qgis_data/ms/MSCities_Geo_Pts.shp", "MSCities_Geo_Pts", "ogr")
    
  2. Next, load the polygon layer:
    lyrPoly = QgsVectorLayer("/qgis_data/ms/GIS_CensusTract_poly.shp", "GIS_CensusTract_poly", "ogr")
    
  3. Add the layers to the map using a list:
    QgsMapLayerRegistry.instance().addMapLayers([lyrPts,lyrPoly])
    
  4. Access the polygon layer's features:
    ftsPoly = lyrPoly.getFeatures()
    
  5. Now, iterate through the polygon's features:
    for feat in ftsPoly:
    
  6. Grab each feature's geometry:
     geomPoly = feat.geometry() 
    
  7. Access the point features and filter the point features by the polygon's bounding box:
     featsPnt = lyrPts.getFeatures(QgsFeatureRequest().setFilterRect(geomPoly.boundingBox()))
    
  8. Iterate through each point and check whether it's within the polygon itself:
     for featPnt in featsPnt:
     if featPnt.geometry().within(geomPoly):
    
  9. If the polygon contains the point, print the point's ID and select the point:
     print featPnt.id()
     lyrPts.select(featPnt.id())
    
  10. Now, set the polygon layer as the active map layer:
    iface.setActiveLayer(lyrPoly)
    
  11. Zoom to the polygon layer's maximum extent:
    iface.zoomToActiveLayer()
    

Verify that your map looks similar to the following image:

How to do it...

How it works...

While QGIS has a number of tools for spatial selection, PyQGIS doesn't have a dedicated API for these types of functions. However, there are just enough methods in the API, thanks to the underlying ogr/GEOS library, that you can easily create your own spatial filters for two layers. Step 7 isn't entirely necessary, but we gain some efficiency using the bounding box of the polygon to limit the number of point features we're examining. Calculations involving rectangles are far quicker than detailed point-in-polygon queries. So, we quickly reduce the number of points we need to iterate through for the more expensive spatial operations.

主站蜘蛛池模板: 青浦区| 德钦县| 清新县| 汉沽区| 咸丰县| 兴安盟| 吕梁市| 芮城县| 永宁县| 巫山县| 当阳市| 泗阳县| 怀柔区| 勐海县| 上栗县| 靖边县| 宁安市| 南木林县| 天全县| 临安市| 巴彦淖尔市| 余干县| 安康市| 焉耆| 天祝| 孝昌县| 奇台县| 宿松县| 年辖:市辖区| 尤溪县| 望城县| 武冈市| 惠东县| 县级市| 凌云县| 介休市| 霞浦县| 余干县| 台北县| 聂荣县| 涡阳县|