How to use GDAL to save raster files in Python

Fellow researchers and open-source GIS enthusiasts,

Welcome to my blog!

I’d like to start with a disclaimer – I may be a researcher of this very area but that doesn’t mean everything I do or write here will work for you, in your own desktop configurations and package versions. I have no responsibility if you lose data or mess up your installation. I also do not authorize any copies of my content.

Today I will write about how to save rasters using GDAL in Python, without using rasterio.

I also have a previous post about reading and saving raster data on Python, that can provide useful context to this one.

Solutions using rasterio work, we all know that. However, installing rasterio and GDAL in the same conda environment is many times a challenge… So I usually try not to install both of them together.

I am showing to you today what I think is, in my experience, the easiest way to save raster data using GDAL.

For reference, I am using Python 3.10, GDAL 3.0, and my IDE of choice is Spyder 5.3.3, running on a Windows 10 machine.

Import GDAL package in Python

from osgeo import gdal

Import OGR Spatial Reference

import osr

Choose to save as .tif

driver = gdal.GetDriverByName('GTiff')

Define the location to save your raster

path_output = 'C:/path_to_save.tif'

Define the number of pixels in your raster

Usually, it is defined by the size of the array you want to save as a raster.

size1,size2=img_as_array.shape

Create the space where your dataset will be placed

Attention: the options for your driver might be different from those.

Specifically, be attentive to the data type of the array you are saving as a raster (the fifth argument). Please check all available data types on the GDAL documentation.

The fourth argument, here simply listed as “1”, is the number of bands your output raster should have.

dataset_output=driver.Create(path_output,size2,size1,1,gdal.GDT_Float64)

It is also possible to append other options, as strings, after these five mandatory arguments.

Learn more about GDAL Create on GDAL docs.

Set the geographic transformation

Usually, this information will be based on another raster.

Use the snippet below if you need to extract a GeoTransform from another raster.

orig_file = gdal.Open(filename, gdal.GA_ReadOnly)
GT_output = orig_file.GetGeoTransform()

Set the GeoTransform:

dataset_output.SetGeoTransform(GT_output)

Set spatial reference for the GeoTIFF

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')

Note: Replace “WGS84” for another spatial reference, when it’s the case.

dataset_output.SetProjection(srs.ExportToWkt())

Finally, write the array into the raster.

export= dataset_output.GetRasterBand(1).WriteArray(img_as_array)

Free the memory the driver was using

export = None

Done!

This should result in a saved raster.

Extras

What if I prefer to use rasterio (a.k.a. rio)?

You can definitely keep using your preferred packages. This post simply shows that it is possible to save rasters on Python using GDAL itself.

Why should we set the variable “export” to None?

I believe GDAL does have an explicit way to close the datasets, GDALClose. But by setting the exported data to None value, we flag that we are no longer using that memory space. That also helps to prevent errors such as trying to open a file in QGIS or ArcMap and receiving an error message that the file is in use by another program, or that the raster could not be opened.

Why should we search the GTiff driver by name?

This is the simplest way to pick the right driver, in order to save your file with the desired filetype.

The list of the available drivers can be found on the GDAL documentation.