Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
    
Try to see this. gis.stackexchange.com/questions/34919/… –  diegogb Apr 13 '13 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 Apr 13 '13 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 Apr 13 '13 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 Apr 13 '13 at 13:59
add comment

4 Answers

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
input.bounds
(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']))})

Result:

enter image description here

share|improve this answer
    
This looks like exactly the procedure i need. I'll have a look at it –  Curlew Apr 13 '13 at 19:44
add comment

ogrinfo -spat:

-spat xmin ymin xmax ymax

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

available in ogr2ogr too.

share|improve this answer
add comment

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

share|improve this answer
    
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 Apr 13 '13 at 14:49
add comment
up vote 0 down vote accepted

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
    '''
    ext=[]
    xarr=[0,cols]
    yarr=[0,rows]

    for px in xarr:
        for py in yarr:
            x=gt[0]+(px*gt[1])+(py*gt[2])
            y=gt[3]+(px*gt[4])+(py*gt[5])
            ext.append([x,y])
        yarr.reverse()
    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
    else:
        # 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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.