How to use the gdal.RasterizeLayer function in GDAL
To help you get started, we’ve selected a few GDAL examples, based on popular ways it is used in public projects.
Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
DigitalGlobe / mltools / mltools / simple_lulc / simple_lulc.py View on Github 
# Create memory target raster
target_ds = gdal.GetDriverByName('MEM').Create('',
xcount, ycount, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((
xmin, pixelWidth, 0,
ymax, 0, pixelHeight,
))
# Create target raster projection
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(raster.GetProjectionRef())
target_ds.SetProjection(raster_srs.ExportToWkt())
# Rasterize zone polygon to raster
gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])
# Read raster as arrays
dataraster = raster.ReadAsArray(xoff, yoff, xcount, ycount)
dataraster = dataraster.astype(np.float)
datamask = target_ds.ReadAsArray(0, 0, xcount, ycount).astype(np.float)
# replicate mask for each band
datamask = np.dstack([datamask for i in range(nbands)])
datamask = datamask.transpose(2,0,1)
# Mask zone of raster
zoneraster = np.ma.masked_array(dataraster, np.logical_not(datamask))
yield (feat, zoneraster)xoff = int((xmin - xOrigin)/pixelWidth)
yoff = int((yOrigin - ymax)/pixelWidth)
xcount = int((xmax - xmin)/pixelWidth)+1
ycount = int((ymax - ymin)/pixelWidth)+1
# Create memory target raster
target_ds = gdal.GetDriverByName('MEM').Create('', xcount, ycount, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((xmin, pixelWidth, 0, ymax, 0, pixelHeight))
# Create for target raster the same projection as for the value raster
raster_srs = osr.SpatialReference();
raster_srs.ImportFromWkt(raster.GetProjectionRef());
target_ds.SetProjection(raster_srs.ExportToWkt());
# Rasterize zone polygon to raster
gdal.RasterizeLayer(target_ds, [1], layer, options = ["ALL_TOUCHED=TRUE", "BURN_VALUE_FROM"]);
# Read raster as arrays
mean = [];
valor = 0
iteration = 0
bands = raster.RasterCount;
#get all rasters
for i in range(1, bands+1):
banddataraster = raster.GetRasterBand(i);
dataraster = banddataraster.ReadAsArray(xoff, yoff, xcount, ycount).astype(numpy.float);
bandmask = target_ds.GetRasterBand(1);
datamask = bandmask.ReadAsArray(0, 0, xcount, ycount).astype(numpy.float);
zoneraster = numpy.ma.masked_array(dataraster, numpy.logical_not(datamask));
for line in zoneraster:
for colum in line:xoff = int((xmin - xOrigin)/pixelWidth)
yoff = int((yOrigin - ymax)/pixelWidth)
xcount = int((xmax - xmin)/pixelWidth)+1
ycount = int((ymax - ymin)/pixelWidth)+1
# Create memory target raster
target_ds = gdal.GetDriverByName('MEM').Create('', xcount, ycount,
1, gdal.GDT_Byte)
target_ds.SetGeoTransform( (xmin, pixelWidth, 0,
ymax, 0, pixelHeight,))
# Create target raster projection
target_ds.SetProjection(targetSR.ExportToWkt())
# Rasterize zone polygon to raster
gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])
# Read raster as arrays
dataraster = raster.ReadAsArray(xoff, yoff, xcount, ycount).astype(np.float)
datamask = target_ds.ReadAsArray(0, 0, xcount, ycount).astype(np.float)
# replicate mask for each band
datamask = np.dstack([datamask for i in range(nbands)])
datamask = datamask.transpose(2,0,1)
# Mask zone of raster
zoneraster = np.ma.masked_array(dataraster, np.logical_not(datamask))
yield ( feat, zoneraster )remove = True
saveto = tempfile.NamedTemporaryFile(suffix='.tif')
saveto = saveto.name
else:
remove = False
# Create the new raster
driver = gdal.GetDriverByName('GTiff')
outRas = driver.Create(saveto, cols, rows, 1)
outRas.SetProjection(lyrSrs)
outRas.SetGeoTransform(geoTrans)
outRas.GetRasterBand(1).Fill(1)
# Rasterize the layer
if method.lower() == 'touches':
gdal.RasterizeLayer(outRas,[1],layer,None, None, [0], ['ALL_TOUCHED=TRUE'])
else: # Just default to this.
gdal.RasterizeLayer(outRas,[1],layer,None, None, [0])
arr = outRas.ReadAsArray()
if remove:
os.remove(saveto)
# Return the numpy array
return arr# Create memory vector layer
mem_ds = ogr.GetDriverByName('Memory').CreateDataSource('out')
mem_layer = mem_ds.CreateLayer(
geom.GetGeometryName(),
None,
geom.GetGeometryType()
)
mem_layer.CreateFeature(feat.Clone())
# Create for target raster the same projection as for the value raster
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(raster.GetProjectionRef())
target_ds.SetProjection(raster_srs.ExportToWkt())
# Rasterize zone polygon to raster
gdal.RasterizeLayer(target_ds, [1], mem_layer, burn_values=[1])
# Read raster as arrays
banddataraster = raster.GetRasterBand(1)
try:
dataraster = banddataraster.ReadAsArray(
xoff, yoff, xcount, ycount).astype(numpy.float)
except AttributeError:
# Nothing within bounds, move on to next polygon
properties = feature['properties']
for p in ['count', 'sum', 'mean', 'median', 'min', 'max', 'stddev']:
properties[p] = None
yield(feature)
else:
# Get no data value of array
noDataValue = banddataraster.GetNoDataValue()
if noDataValue:xcount = int((xmax - xmin)/pixelWidth)+1
ycount = int((ymax - ymin)/pixelWidth)+1
# Create memory target raster
target_ds = gdal.GetDriverByName('MEM').Create('', xcount, ycount, gdal.GDT_Byte)
target_ds.SetGeoTransform((
xmin, pixelWidth, 0,
ymax, 0, pixelHeight,
))
# Create for target raster the same projection as for the value raster
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(raster.GetProjectionRef())
target_ds.SetProjection(raster_srs.ExportToWkt())
# Rasterize zone polygon to raster
gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])
# Read raster as arrays
banddataraster = raster.GetRasterBand(1)
# error here
dataraster = banddataraster.ReadAsArray(xoff, yoff, xcount, ycount).astype(numpy.float)
bandmask = target_ds.GetRasterBand(1)
datamask = bandmask.ReadAsArray(0, 0, xcount, ycount).astype(numpy.float)
# Mask zone of raster
zoneraster = numpy.ma.masked_array(dataraster, numpy.logical_not(datamask))
# Calculate statistics of zonal raster
value = numpy.mean(zoneraster)
print value
return value# Output
out_driver = gdal.GetDriverByName('GTiff')
if os.path.exists(output_tiff):
out_driver.Delete(output_tiff)
out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
1, gdal.GDT_Int16)
out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
out_source.SetProjection(inp_srs.ExportToWkt())
out_lyr = out_source.GetRasterBand(1)
out_lyr.SetNoDataValue(NoData_value)
# Rasterize
if field_name:
gdal.RasterizeLayer(out_source, [1], inp_lyr,
options=["ATTRIBUTE={0}".format(field_name)])
else:
gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])
# Save and/or close the data sources
inp_source = None
out_source = None
# Return
return output_tiffmem_raster = mem_drv.Create(
'',
model_dataset.RasterXSize,
model_dataset.RasterYSize,
1,
dtype
)
mem_raster.SetProjection(model_dataset.GetProjection())
mem_raster.SetGeoTransform(model_dataset.GetGeoTransform())
mem_band = mem_raster.GetRasterBand(1)
mem_band.Fill(nodata_val)
mem_band.SetNoDataValue(nodata_val)
# http://gdal.org/gdal__alg_8h.html#adfe5e5d287d6c184aab03acbfa567cb1
# http://gis.stackexchange.com/questions/31568/gdal-rasterizelayer-doesnt-burn-all-polygons-to-raster
err = gdal.RasterizeLayer(
mem_raster,
[1],
shape_layer,
None,
None,
[1],
options
)
assert err == gdal.CE_None
return mem_raster.ReadAsArray()# Create memory target raster
target_ds = gdal.GetDriverByName('MEM').Create(
'', xcount, ycount, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((
xmin, pixelWidth, 0,
ymax, 0, pixelHeight,
))
# Create for target raster the same projection as for the value raster
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(raster.GetProjectionRef())
target_ds.SetProjection(raster_srs.ExportToWkt())
# Rasterize zone polygon to raster
gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])
# Read raster as arrays
banddataraster = raster.GetRasterBand(1)
try:
dataraster = banddataraster.ReadAsArray(
xoff, yoff, xcount, ycount).astype(numpy.float)
except AttributeError:
# Nothing within bounds, move on to next polygon
properties = feature['properties']
for p in ['count', 'sum', 'mean', 'median', 'min', 'max', 'stddev']:
properties[p] = None
yield(feature)
else:
# Get no data value of array
noDataValue = banddataraster.GetNoDataValue()
if noDataValue:
Popular GDAL functions
- gdal.AllRegister
- gdal.ApplyGeoTransform
- gdal.AutoCreateWarpedVRT
- gdal.GA_ReadOnly
- gdal.GA_Update
- gdal.GCP
- gdal.GCPsToGeoTransform
- gdal.GDT_Byte
- gdal.GDT_Float32
- gdal.GeneralCmdLineProcessor
- gdal.GetDataTypeName
- gdal.GetDriverByName
- gdal.GetLastErrorMsg
- gdal.InvGeoTransform
- gdal.Open
- gdal.PushErrorHandler
- gdal.RasterizeLayer
- gdal.ReprojectImage
- gdal.SetConfigOption
- gdal.UseExceptions