PyQGIS: Create Raster

In a previous tutorial I showed you how to access raster values and data with the Geospatial Data Abstraction Library (GDAL). This tutorial will show you how to create a raster with GDAL.

Start by importing the gdal and osr Python modules. The osr module is used for handling spatial references. Also import numpy.

from osgeo import gdal, osr
import numpy as np

Now define a variable containing the path to the output file you will create.

fn = 'path/to/raster.tif'

To create a raster we need to have data to write to the file. Here we’ll create a data that represents a 2-dimensional raster with 10 rows and 10 columns. The 2D array will be filled with zeros.

rasterband = np.zeros((10,10))

GDAL uses drivers to create new files. There are different drivers for different file formats. Generally, I use geotiff, or tif, files to store my raster data. GDAL supports a number of other drivers. The code below retrieves a geotiff driver then uses the driver to create a new raster dataset, ds, which has 10 rows, 10 columns, 1 band, and a 32 bit float data type.

driver = gdal.GetDriverByName('GTiff')
ds = driver.Create(fn, xsize=10, ysize=10, bands=1, eType=gdal.GDT_Float32)

You’ve now created a raster dataset, but it is empty and does not contain any data. You can add data be retrieving a raster band and writing data to that band. The WriteArray() method allows us to write to a raster band with a numpy array that matches the raster’s dimensions.

ds.GetRasterBand(1).WriteArray(rasterband)

The raster dataset now contains data, but has no reference to a geographic location. To spatially reference a raster you need to provide two pieces of information: the affine geotransformation, which gives the top left coordinates and defines the cell size, and the spatial reference, which defines the projection.

With GDAL the geotransformation is defined as a list with the following values [top-left x coord, cell width, 0, top-left y coord, 0, cell height]. Note that cell height is negative for north-up rasters (nearly all rasters are north-up). The zero values are for rotation and shearing values. Generally, these are set to zero. The geotransform information below is for a UTM projection, so the units for coordinates and cell dimensions are in meters. Once the geotransform is defined it is set by calling SetGeoTransform().

geot = [500000, 10, 0, 4600000, 0, -10]
ds.SetGeoTransform(geot)

Now the spatial reference (also referred to as coordinate reference system or projection) can be defined for the raster data set. This is where you’ll use the osr module. First create a spatial reference object. For the Universal Transverse Mercator the SetUTM(Zone, North=1 or South=2)method is called. Other methods can set the spatial reference from well-known text or EPSG code. Then the Geographic Coordinate System is set. In this case the North American Datum of 1983. Then we can set the projection for our raster dataset. The SetProjection() method requires the spatial reference to be input in well-known text format.

srs = osr.SpatialReference()
srs.SetUTM(12,1)
srs.SetWellKnownGeogCS('NAD83')
ds.SetProjection(srs.ExportToWkt())

The last thing to do is to close the rater datasource by setting it equal to None. Then you can add the new raster to the QGIS interface.

ds = None
rlayer = iface.addRasterLayer(fn)

Below is the full code block we developed. You can also watch the video at the end of this post to see a step-by-step demonstration.

from osgeo import gdal, osr
import numpy as np

fn = 'path/to/raster.tif'

rasterband = np.zeros((10,10))

driver = gdal.GetDriverByName('GTiff')
ds = driver.Create(fn, xsize=10, ysize=10, bands=1, eType=gdal.GDT_Float32)
ds.GetRasterBand(1).WriteArray(rasterband)

geot = [500000, 10, 0, 4600000, 0, -10]
ds.SetGeoTransform(geot)
srs = osr.SpatialReference()
srs.SetUTM(12,1)
srs.SetWellKnownGeogCS('NAD83')
ds.SetProjection(srs.ExportToWkt())
ds = None

rlayer = iface.addRasterLayer(fn)