Exporting data to GeoTIFF files

  • Compatability: Notebook currently compatible with both the NCI and DEA Sandbox environments

  • Products used: ga_ls8c_ard_3


At the end of an analysis, it can be useful to export data to a GeoTIFF (e.g. outputname.tif) file, either to save results or to allow for exploring results in a GIS software platform (e.g. ArcGIS or QGIS).


This notebook shows a number of ways to export a GeoTIFF file:

  1. Exporting a single-band, single time-slice GeoTIFF from an xarray object loaded through a dc.load query

  2. Exporting a multi-band, single time-slice GeoTIFF from an xarray object loaded through a dc.load query

  3. Exporting multiple GeoTIFFs, one for each time-slice of an xarray object loaded through a dc.load query

  4. Exporting a single-band GeoTIFF when handling a simple array and/or an xarray that doesn’t contain transform, crs or affine metadata

Getting started

To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.

Load packages

%matplotlib inline

import datacube
from datacube.helpers import write_geotiff
import numpy as np
import xarray as xr
import sys

import dea_datahandling
from dea_plotting import rgb

Connect to the datacube

dc = datacube.Datacube(app='Exporting_GeoTIFFs')

Load Landsat 8 data from the datacube

Here we are loading in a timeseries of Landsat 8 satellite images through the datacube API. This will provide us with some data to work with.

# Create a query object
query = {
    'x': (153.40, 153.45),
    'y': (-28.85, -28.90),
    'time': ('2018-01', '2018-03'),
    'measurements': ['nbart_red', 'nbart_green', 'nbart_blue'],
    'output_crs': 'EPSG:3577',
    'resolution': (-30, 30),
    'group_by': 'solar_day'

# Load available data from the Landsat 8
ds = dc.load(product='ga_ls8c_ard_3', **query)

# Print output data

Dimensions:      (time: 6, x: 191, y: 212)
  * time         (time) datetime64[ns] 2018-01-03T23:42:39.393181 ... 2018-03-24T23:42:01.382412
  * y            (y) float64 -3.307e+06 -3.307e+06 ... -3.313e+06 -3.313e+06
  * x            (x) float64 2.053e+06 2.053e+06 ... 2.058e+06 2.058e+06
Data variables:
    nbart_red    (time, y, x) int16 411 422 404 396 509 ... 244 285 307 251 247
    nbart_green  (time, y, x) int16 625 640 595 580 697 ... 343 407 470 372 371
    nbart_blue   (time, y, x) int16 435 430 421 412 542 ... 273 308 320 274 270
    crs:      EPSG:3577

Plot an rgb image to confirm we have data

The white regions are cloud cover.

rgb(ds, index=5, percentile_stretch=(0.1, 0.8))


Export a single-band, single time-slice GeoTIFF

This method uses the datacube.helpers function write_geotiff to export a simple single-band, single time-slice GeoTIFF.

Note: An important caveat to using this function is that it requires an xarray.Dataset, supplying an xarray.DataArray will return an error.

# Select a single time-slice and a single band from the dataset.
# The [['nbart_red']] keeps the array as a `xarray.Dataset`.
singleBandtiff = ds[['nbart_red']].isel(time=1)

# Write GeoTIFF to a location
write_geotiff('red_band.tif', singleBandtiff)

Export a multi-band, single time-slice GeoTIFF

Here we select a single time and export all the bands in the dataset using the datacube.helpers.write_geotiff function.

# Select a single time-slice
rgb_tiff = ds.isel(time=1)

# Write multi-band GeoTIFF to a location
write_geotiff('rgb.tif', rgb_tiff)

Export multiple GeoTIFF, one for each time-slice of an xarray

If we want to export all of the time steps in a dataset as a GeoTIFF, we can wrap our write_geotiff function in a for-loop.

for i in range(len(ds.time)):

    # We will use the date of the satellite image to name the GeoTIFF
    date = str(ds.isel(time=i).time.data)[:-19]
    print(f'Writing {date}')

    # Write GeoTIFF
    write_geotiff(f'{date}.tif', ds.isel(time=i))

Writing 2018-01-03
Writing 2018-01-19
Writing 2018-02-04
Writing 2018-02-20
Writing 2018-03-08
Writing 2018-03-24

Exporting simple arrays without geospatial metadata

Often in a workflow we create a numpy array that doesn’t have the metadata required to work with functions like datacube.helpers.write_geotiff. In these instances, we can use the function dea_datahandling.array_to_geotiff to export numpy arrays that don’t have metadata attached.

First we will perform some trivial analysis that will require shifting our datacube-acquired xarray object into a numpy array, then we will export our numpy array as a GeoTIFF. This approach will work on any numpy array as long as we know the projection and the coordinates (which we can get from the original xarray.Dataset).

# Convert one of the xarray.DataArrays inside the xarray.Dataset into a numpy array
arr = ds.nbart_red.isel(time=1).values

# Perform some trivial analysis
simple_array = np.where(arr >= 1, np.nan, 5)

Export a single band GeoTIFF

We will use the function dea_datahandling.array_to_geotiff to export our simple array:

# First grab the spatial information from our datacube xarray object
transform = ds.geobox.transform.to_gdal()
projection = ds.geobox.crs.wkt

# Export the array

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube Slack channel or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on Github.

Last modified: December 2019

Compatible datacube version:



Browse all available tags on the DEA User Guide’s Tags Index

Tags: sandbox compatible, NCI compatible, GeoTIFF, dc.load, landsat 8, array_to_geotiff, datacube.helpers.write_geotiff, exporting data