2

I have a raster layer (coloured squares) and a vector layer that represents water bodies (blue coloured polygones; see image below). I am looking for a way to exclude / remove the cells that are covered by more than a certain percentage of water (say 60%). How can I do it using gdal and / or python?

enter image description here

1 Answer 1

2

Next code can be used to do that:

from osgeo import gdal, osr

registry = QgsMapLayerRegistry.instance()

raster = registry.mapLayersByName('my_raster')
provider = raster[0].dataProvider()
filename = provider.dataSourceUri()
water_bodies = registry.mapLayersByName('water_bodies_mult')

dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)
transform = dataset.GetGeoTransform()
data = band.ReadAsArray()

feat_wb = water_bodies[0].getFeatures().next() 

extent = raster[0].extent()

xmin, ymin, xmax, ymax = extent.toRectF().getCoords()

xsize = raster[0].rasterUnitsPerPixelX()
ysize = raster[0].rasterUnitsPerPixelY()

pixel_area = xsize*ysize

rows = raster[0].height()
cols = raster[0].width()

x1 = xmin
y1 = ymax

for i in range(rows):
    for j in range(cols):
        p1 = QgsPoint(x1, y1)
        p2 = QgsPoint(x1 + xsize, y1 - ysize)
        rect = QgsRectangle(p1, p2)
        geom = QgsGeometry.fromWkt(rect.asWktPolygon())
        if geom.intersects(feat_wb.geometry()):
            new_geom = geom.intersection(feat_wb.geometry())
            perc = (new_geom.area()/pixel_area)*100
            if perc >= 60.:
                data[i][j] = -9999

        x1 += xsize
    x1 = xmin
    y1 -= ysize

# Create gtif file 
driver = gdal.GetDriverByName("GTiff")

output_file = "/home/zeito/pyqgis_data/new_raster.tif"

dst_ds = driver.Create(output_file, 
                       cols, 
                       rows, 
                       1, 
                       gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( data )

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(transform)

wkt = dataset.GetProjection()

# setting spatial reference of output raster 
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
dst_ds.GetRasterBand(1).SetNoDataValue(-9999)

#Close output raster dataset 
dataset = None
dst_ds = None

For a situation represented by next image:

enter image description here

after running above code and load resulting raster (new_raster.tif), it can be observed, at next image, that the cells that were covered by more than 60 % of water were excluded (resulting raster is transparent in these areas).

enter image description here

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.