Friday 1 June 2018

Determining if shapefile and raster overlap in Python using OGR/GDAL?



I am building a script in python using OGR/GDAL.


I have a set of shapefiles and a set of GeoTiff raster files.


I would like to have my script ignore shapefiles if they do not intersect with the raster area.


The shapefile is not a rectangle, so I cannot simply compare the xmin/xmax,ymin/ymax values returned by layer.GetExtent(). I need the actual polygon representing it's overall shape, and then some way of determining if that polygon intersects with the raster square.


I was thinking I could somehow merge all the polygons in the shapefile into one feature, and then read the geometry on that feature, and then compare that information to the raster extent. However, I am unsure of specifically how to execute this.



  1. How to extract border polygon information from shapefile?

  2. How to determine if that polygon intersects a given square area?




Answer



The following script determines the bounding box of a raster and creates based on the bounding box a geometry.


import ogr, gdal

raster = gdal.Open('sample.tif')
vector = ogr.Open('sample.shp')

# Get raster geometry
transform = raster.GetGeoTransform()
pixelWidth = transform[1]

pixelHeight = transform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize

xLeft = transform[0]
yTop = transform[3]
xRight = xLeft+cols*pixelWidth
yBottom = yTop-rows*pixelHeight

ring = ogr.Geometry(ogr.wkbLinearRing)

ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xLeft, yTop)
rasterGeometry = ogr.Geometry(ogr.wkbPolygon)
rasterGeometry.AddGeometry(ring)

Next, the geometry of the vector polygon is determined. This answers your first question.


# Get vector geometry

layer = vector.GetLayer()
feature = layer.GetFeature(0)
vectorGeometry = feature.GetGeometryRef()

Last, the geometry of the vector and raster are tested for intersection (returns True or False). This answers your second question.


print rasterGeometry.Intersect(vectorGeometry)

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...