1

I need documentation on how to apply a spatial filter (using OGR's SetSpatialFilter).

Setting the feature layer by which another layer should be filtered seems rather straightforward. But, how do you set the type of spatial filtering ('intersect', 'contain' etc.). And where do you find a list of optional types?

In the case at hand, I would like to filter out features from one layer (let's call it 'TilePolys')that are touched ('intersects with' (?)) by features from another layer (let's call it 'Buildings).

2

1 Answer 1

4

Here is a link to the documentation, with the following content:

SetSpatialFilter(Layer self, int iGeomField, Geometry filter)

void OGR_L_SetSpatialFilter(OGRLayerH hLayer, OGRGeometryH hGeom)

Set a new spatial filter.

This function set the geometry to be used as a spatial filter when
fetching features via the OGR_L_GetNextFeature() function. Only
features that geometrically intersect the filter geometry will be
returned.

Currently this test is may be inaccurately implemented, but it is
guaranteed that all features whose envelope (as returned by
OGR_G_GetEnvelope()) overlaps the envelope of the spatial filter will
be returned. This can result in more shapes being returned that should
strictly be the case.

Starting with GDAL 2.3, features with null or empty geometries will
never be considered as matching a spatial filter.

This function makes an internal copy of the passed geometry. The
passed geometry remains the responsibility of the caller, and may be
safely destroyed.

For the time being the passed filter geometry should be in the same
SRS as the layer (as returned by OGR_L_GetSpatialRef()). In the future
this may be generalized.

This function is the same as the C++ method
OGRLayer::SetSpatialFilter.

Parameters:
-----------

hLayer:  handle to the layer on which to set the spatial filter.

hGeom:  handle to the geometry to use as a filtering region. NULL may
be passed indicating that the current spatial filter should be
cleared, but no new one instituted.  ```


Using python would be as simple as follows:

from osgeo import ogr

ds = ogr.Open(shapefile)
lyr = ds.GetLayer()  # TilePolys

wkt = "POLYGON ((-103.81402655265633 50.253951270672125,-102.94583419409656 51.535568561879401,-100.34125711841725 51.328856095555651,-100.34125711841725 51.328856095555651,-93.437060743203844 50.460663736995883,-93.767800689321859 46.450441890315041,-94.635993047881612 41.613370178339181,-100.75468205106476 41.365315218750681,-106.12920617548238 42.564247523428456,-105.96383620242338 47.277291755610058,-103.81402655265633 50.253951270672125))"

geometry = ogr.CreateGeometryFromWkt(wkt)   # Buildings

lyr.SetSpatialFilter(geometry)

In this case the geometry is made from a Wkt but you can use whatever geometry you like. Usually you would import the shapefile where your 'filtering' feature(s) is (are) and use it to set the spatial filter.

As the comment above says, intersection is the only filter that works with SetSpatialFilter and would work in your case.

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.