Buffer

CreateBuffer.py

import os
import gdal, ogr, osr
import numpy as np

def createBuffer(inputfn, outputBufferfn, bufferDist, srs):

    inputds = ogr.Open(inputfn)
    inputlyr = inputds.GetLayer()

    shpdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputBufferfn):
        shpdriver.DeleteDataSource(outputBufferfn)
    outputBufferds = shpdriver.CreateDataSource(outputBufferfn)
    bufferlyr = outputBufferds.CreateLayer(outputBufferfn,
                                           srs = srs,
                                           geom_type=ogr.wkbPolygon)
    featureDefn = bufferlyr.GetLayerDefn()

    for feature in inputlyr:
        ingeom = feature.GetGeometryRef()
        geomBuffer = ingeom.Buffer(bufferDist)

        outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(geomBuffer)
        bufferlyr.CreateFeature(outFeature)

# Read input data
# NI LC map
d = gdal.Open( 'lcm2007_25m_ni.tif' )
b = d.GetRasterBand(1)
proj = d.GetProjection()
gt = d.GetGeoTransform()

# Water
LC_ID = 16

mask = b.ReadAsArray()
mask = np.where( mask  == LC_ID, 1, 0)
rows, cols = mask.shape

# Create MEM dataset
driver = gdal.GetDriverByName('MEM')
mask_ds = driver.Create('', cols, rows, 1)
mask_band = mask_ds.GetRasterBand(1)
mask_band.WriteArray( mask )

# Get projection information
srs = osr.SpatialReference()
srs.ImportFromWkt( proj )

# Create tmp.shp
dst_layername = "tmp_poly.shp"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername )
dst_layer = dst_ds.CreateLayer( dst_layername, srs = srs )
newField = ogr.FieldDefn( 'ID', ogr.OFTInteger )
dst_layer.CreateField( newField )

# Vectorize
LC_ID = 7
gdal.Polygonize( b, mask_band, dst_layer, 0, [], callback=None )
dst_ds = None

# Create buffer
dst_buffer_layername = "buffer.shp"
bufferDist = 100.0 # in meters
createBuffer(dst_layername, dst_buffer_layername, bufferDist, srs)

# Rasterize buffer
buffer_vector = ogr.Open( dst_buffer_layername )
buffer_layer = buffer_vector.GetLayer()

output = 'buffer.tif'
driver = gdal.GetDriverByName('GTiff')
buffer_ds = driver.Create( output, cols, rows, 1, gdal.GDT_Byte)
buffer_ds.SetProjection( proj )
buffer_ds.SetGeoTransform( gt )
band = buffer_ds.GetRasterBand(1)
NoData_value = 0
band.SetNoDataValue( NoData_value )
band.FlushCache()
gdal.RasterizeLayer( buffer_ds, [1], buffer_layer )