# Function to create a new 1-band raster and use projection
# and geotransform from another raster.
from osgeo import gdal
def make_raster(in_ds, fn, data, data_type, nodata=None):
"""Create a one-band GeoTIFF.
in_ds - datasource to copy projection and geotransform from
fn - path to the file to create
data - NumPy array containing data to write
data_type - output data type
nodata - optional NoData value
"""
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(
fn, in_ds.RasterXSize, in_ds.RasterYSize, 1, data_type)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)
if nodata is not None:
out_band.SetNoDataValue(nodata)
out_band.WriteArray(data)
out_band.FlushCache()
out_band.ComputeStatistics(False)
return out_ds
