Friday, 10 February 2017

python - Get vector features inside a specific extent


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



Answer



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



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...