The problem:

I have a very large vector layer with many features and a much smaller raster layer in a defined region. I want to get only those vector features which are within the raster layers extent (extent = (xmin,xmax,ymin,ymax) ).

Is there anything like a standard SQL-query (something like SELECT * FROM layer WHERE EXTENT < extent) or another command (some test?) i could use to get only those features that are within a given extent?

EDIT: Added python code to do a bounding box intersection test for vector features using just gdal and ogr python binding

  • Try to see this. gis.stackexchange.com/questions/34919/…
    – diegogb
    Commented Apr 13, 2013 at 10:09
  • Sorry, but i don't use POSTGIS and as far as i know the ST_Within parameter works only with postgis and not with ogr sql queries
    – Curlew
    Commented Apr 13, 2013 at 10:20
  • That is essentially what RFC-6 (trac.osgeo.org/gdal/wiki/rfc6_sqlgeom) provided, and you've linked to. You have the operators (<, > and BETWEEN), so what is the problem?
    – BradHards
    Commented Apr 13, 2013 at 11:25
  • ahh, i still fail to get it. Let's say i have a tuple containing an extent of a layer and i want all vector features (partly)within this extent in python. How would you do this? Can you give me some pseudo-code?
    – Curlew
    Commented Apr 13, 2013 at 13:59

5 Answers 5


You don't need a SQL-query to do that, only Python with, once again, the modules Fiona and Shapely of Sean Gillies.

I want only the records which are within the blue frame (analogy of a raster layer).

enter image description here

See the Fiona user Manual, the filter() method returns an iterator over records that intersect a given (minx, miny, maxx, maxy) bounding box:

from shapely.geometry import mapping, shape
import fiona
# Read the original Shapefile
input = fiona.open('data.shp', 'r')
# bounds of the original shapefile
(258018.9133083854, 158162.863836, 268763.670357, 162621.686305)
# clip the shapefile with the raster bounds 
clipped = input.filter(bbox=((262236.3101588468, 159973.80344954136, 263491.7250217228, 160827.485556297)))
# create the clipped shapefile with the same schema
clipped_schema = input.schema.copy()
with fiona.collection('clipped.shp', 'w', 'ESRI Shapefile', clipped_schema) as output:
    for elem in clipped:
           output.write({'properties': elem['properties'],'geometry': mapping(shape(elem['geometry']))})


enter image description here

  • This looks like exactly the procedure i need. I'll have a look at it
    – Curlew
    Commented Apr 13, 2013 at 19:44

For those who want the box to have hard, cropped borders (rather than selecting features overlapping or falling within the box), I found a very easy way to do this with GeoPandas. I use an unzipped version of the countries file available here.

import geopandas as gp
from shapely.geometry import Polygon

# Point to the folder containing the unzipped shapefile
basemap = gp.read_file('ne_10m_admin_0_countries')

# Set the bounds of the crop box
xmin, xmax, ymin, ymax = 6, 18, 36, 50 
bounds = Polygon( [(xmin,ymin), (xmin, ymax), (xmax, ymax), (xmax,ymin)] )

# Crop all polygons and take the part inside the bounding box
basemap['geometry'] = basemap['geometry'].intersection(bounds)

# Export non-empty geometries to shp
basemap[basemap.geometry.area>0].to_file("countries_CROPPED.shp", driver='ESRI Shapefile')

I extracted Italy from here:
enter image description here

To here:
enter image description here


ogrinfo -spat:

-spat xmin ymin xmax ymax

see: http://gdal.org/ogrinfo.html

available in ogr2ogr too.


here is my current python solution. First i extract the rasters extent and then i loop through the extents of every vector feature and do a bounding box intersection test. I favor this solution because i don't have to install any more dependencies than the standard gdal/ogr bindings.

# Code from here http://gis.stackexchange.com/questions/57834/how-to-get-raster-corner-coordinates-using-python-gdal-bindings
def GetExtent(self,gt,cols,rows):
    ''' Return list of corner coordinates from a geotransform

        @type gt:   C{tuple/list}
        @param gt: geotransform
        @type cols:   C{int}
        @param cols: number of columns in the dataset
        @type rows:   C{int}
        @param rows: number of rows in the dataset
        @rtype:    C{[float,...,float]}
        @return:   coordinates of each corner

    for px in xarr:
        for py in yarr:
    return ext

# Bounding Box intersection test
def BBoxIntersect(self,rasE,polyE):
    # Get upper left point + height and width of both bbs
    b1_x = rasE[0]
    b1_y = rasE[2]
    b1_w = rasE[1] - rasE[0]
    b1_h = rasE[0] - rasE[2]

    b2_x = polyE[0]
    b2_y = polyE[2]
    b2_w = polyE[1] - polyE[0]
    b2_h = polyE[0] - polyE[2]

    # is b1 on the right side of b2? # is b1 under b2? # is b2 on the right side of b1? # is b2 under b1?
    if (b1_x > b2_x + b2_w- 1 ) or (b1_y > b2_y + b2_h - 1 ) or (b2_x > b1_x + b1_w - 1 ) or (b2_y > b1_y + b1_h - 1 ):                
        # no collision
        return False
        # collision
        return True

import ogr, gdal
# Then do this:
srcImage = gdal.Open("path to your raster")
band = srcImage.GetRasterBand(1) # Get band 1
geoTrans = srcImage.GetGeoTransform() # Get geotransform information
# Now extract the rasters extent
ext = GetExtent(geoTrans,srcImage.RasterXSize,srcImage.RasterYSize)

# Then open your layer and loop through the features.
# Conduct a boundingbox intersection test for every feature
shapef = ogr.Open("path to your vector")
lyr = shapef.GetLayer()

for i in xrange(0,lyr.GetFeatureCount()):
    poly = lyr.GetFeature(i)
    geom = poly.GetGeometryRef()
    f_coord = geom.GetEnvelope() # feature bounding box
    ints = BBoxIntersect(ext,f_coord) # Is the features BB intersecting the raster BB ?
    if ints:
      print poly.getFID(), "is intersecting with the rasters extent"

This works at least for plain boundings boxes, but i admit that depending on the shape of the vector features (non rectangular features) this could result in overlays where there are no raster values beneath.
For instance see this screenshot for problematic vector features. The red colored areas are inside of the features bounding box and could contain raster-values. Therefore the intersection test will result in a TRUE -> Intersection, although the whole feature is outside of the raster. In this case stick with the fiona solution posted by Gene here.

enter image description here


Well, if you need to execute only once, an easy way to clip the vector layer is get the extent of the raster with gdalinfo and then use the -clipsrc [xmin ymin xmax ymax] parameter of ogr2ogr to make the clip

  • Thanks, but i want something like a test (or a subset of features) in python if a vector feature is inside a a given extent. It is gonna be executed multiple times and is used in a scripted context.
    – Curlew
    Commented Apr 13, 2013 at 14:49

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.