Exporting cloud-optimised GeoTIFF files

• Sign up to the DEA Sandbox to run this notebook interactively from a browser

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

• Products used: ga_ls8c_ard_3

Background

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

A Cloud Optimized GeoTIFF (COG) is a regular GeoTIFF file (i.e. that can be opened by GIS software like QGIS or ArcMap) aimed at being hosted on a HTTP file server, with an internal organization that enables more efficient workflows on the cloud.

Description

This notebook shows a number of ways to export a GeoTIFF file using the datacube.utils.cog function write_cog:

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

In addition, the notebook demonstrates several more advanced applications of write_cog: 1. Exporting data from lazily loaded dask arrays 2. Passing in custom rasterio parameters to override the function’s defaults

---

Getting started

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

Load packages

import rasterio
import datacube
from datacube.utils.cog import write_cog

import sys
sys.path.insert(1, '../Tools/')
from dea_tools.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': (149.06, 149.17),
    'y': (-35.27, -35.32),
    'time': ('2020-01', '2020-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
ds

<xarray.Dataset>
Dimensions:      (time: 11, y: 229, x: 356)
Coordinates:
  * time         (time) datetime64[ns] 2020-01-07T23:56:40.592891 ... 2020-03...
  * y            (y) float64 -3.955e+06 -3.955e+06 ... -3.962e+06 -3.962e+06
  * x            (x) float64 1.544e+06 1.544e+06 ... 1.554e+06 1.554e+06
    spatial_ref  int32 3577
Data variables:
    nbart_red    (time, y, x) int16 1691 1729 1738 1734 1695 ... 799 799 819 753
    nbart_green  (time, y, x) int16 1380 1386 1405 1398 1377 ... 879 911 877 808
    nbart_blue   (time, y, x) int16 1170 1168 1175 1171 1156 ... 670 670 685 642
Attributes:
    crs:           EPSG:3577
    grid_mapping:  spatial_ref

xarray.Dataset

• Dimensions:
  • time: 11
  • y: 229
  • x: 356
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2020-01-07T23:56:40.592891 ... 2...

    units :

    seconds since 1970-01-01 00:00:00

    array(['2020-01-07T23:56:40.592891000', '2020-01-16T23:50:27.927367000',
           '2020-01-23T23:56:36.765939000', '2020-02-01T23:50:22.661832000',
           '2020-02-08T23:56:31.774346000', '2020-02-17T23:50:19.161045000',
           '2020-02-24T23:56:28.025939000', '2020-03-04T23:50:13.948244000',
           '2020-03-11T23:56:21.709961000', '2020-03-20T23:50:06.289323000',
           '2020-03-27T23:56:12.950767000'], dtype='datetime64[ns]')

  • y
    (y)
    float64
    -3.955e+06 ... -3.962e+06

    units :

    metre

    resolution :

    -30.0

    crs :

    EPSG:3577

    array([-3955425., -3955455., -3955485., ..., -3962205., -3962235., -3962265.])

  • x
    (x)
    float64
    1.544e+06 1.544e+06 ... 1.554e+06

    units :

    metre

    resolution :

    30.0

    crs :

    EPSG:3577

    array([1543665., 1543695., 1543725., ..., 1554255., 1554285., 1554315.])

  • spatial_ref
    ()
    int32
    3577

    spatial_ref :

    PROJCS["GDA94 / Australian Albers",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",132],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3577"]]

    grid_mapping_name :

    albers_conical_equal_area

    array(3577, dtype=int32)

• Data variables: (3)
  • nbart_red
    (time, y, x)
    int16
    1691 1729 1738 1734 ... 799 819 753

    units :

    1

    nodata :

    -999

    crs :

    EPSG:3577

    grid_mapping :

    spatial_ref

    array([[[1691, 1729, 1738, ..., 1180, 1347, 1175],
            [1680, 1731, 1763, ..., 1140, 1220, 1158],
            [1737, 1778, 1816, ..., 1070, 1047, 1131],
            ...,
            [1766, 1808, 1830, ..., 1175, 1069, 1050],
            [1748, 1731, 1687, ..., 1423, 1193, 1119],
            [1520, 1802, 1868, ..., 1417, 1340, 1296]],
    
           [[4212, 4279, 4452, ..., 5406, 5616, 5215],
            [3856, 4015, 4171, ..., 5230, 5361, 5155],
            [4199, 4290, 4576, ..., 4772, 4718, 4711],
            ...,
            [2714, 2623, 2422, ..., 5703, 5535, 5509],
            [2639, 2544, 2489, ..., 5698, 5615, 5508],
            [2438, 2534, 2468, ..., 5610, 5597, 5633]],
    
           [[1362, 1371, 1395, ...,  969, 1185,  939],
            [1332, 1367, 1392, ...,  939, 1053,  919],
            [1347, 1410, 1462, ...,  813,  749,  851],
            ...,
    ...
            ...,
            [1226, 1705, 2317, ..., 3054, 3378, 3864],
            [1399, 1936, 2214, ..., 3412, 3548, 3919],
            [1399, 1862, 2294, ..., 4103, 4105, 3997]],
    
           [[ 662,  492,  508, ...,  718,  875,  948],
            [ 585,  505,  492, ...,  674,  710,  729],
            [ 453,  445,  461, ...,  564,  494,  567],
            ...,
            [1629, 1840, 1572, ...,  389,  324,  269],
            [1535, 1496, 1709, ...,  449,  474,  369],
            [1370, 1919, 1923, ...,  513,  554,  564]],
    
           [[1161,  868,  689, ...,  857,  957, 1001],
            [1005,  749,  621, ...,  886,  898,  976],
            [ 758,  602,  535, ...,  910,  837,  964],
            ...,
            [ 987,  981,  933, ...,  662,  556,  566],
            [1186,  944,  911, ...,  649,  716,  702],
            [2028, 1222, 1024, ...,  799,  819,  753]]], dtype=int16)

  • nbart_green
    (time, y, x)
    int16
    1380 1386 1405 1398 ... 911 877 808

    units :

    1

    nodata :

    -999

    crs :

    EPSG:3577

    grid_mapping :

    spatial_ref

    array([[[1380, 1386, 1405, ..., 1068, 1162, 1070],
            [1377, 1408, 1419, ..., 1053, 1100, 1079],
            [1397, 1423, 1445, ...,  993, 1000, 1052],
            ...,
            [1630, 1655, 1635, ..., 1174, 1131, 1122],
            [1628, 1543, 1540, ..., 1245, 1187, 1161],
            [1464, 1644, 1636, ..., 1284, 1256, 1242]],
    
           [[4156, 4198, 4302, ..., 5346, 5500, 5163],
            [3671, 3779, 3995, ..., 5145, 5262, 5053],
            [4065, 4115, 4396, ..., 4703, 4614, 4625],
            ...,
            [2717, 2651, 2443, ..., 5618, 5482, 5433],
            [2662, 2525, 2502, ..., 5584, 5493, 5408],
            [2474, 2530, 2478, ..., 5517, 5483, 5508]],
    
           [[1013, 1015, 1031, ...,  764,  911,  754],
            [ 991, 1010, 1029, ...,  730,  788,  714],
            [1000, 1040, 1069, ...,  662,  621,  687],
            ...,
    ...
            ...,
            [1249, 1648, 2240, ..., 2996, 3328, 3785],
            [1379, 1937, 2120, ..., 3344, 3475, 3783],
            [1363, 1848, 2283, ..., 3993, 3999, 3893]],
    
           [[ 784,  751,  769, ...,  605,  754,  734],
            [ 721,  753,  721, ...,  592,  610,  600],
            [ 735,  745,  773, ...,  517,  447,  493],
            ...,
            [1577, 1575, 1428, ...,  475,  344,  262],
            [1461, 1327, 1377, ...,  553,  570,  412],
            [1137, 1568, 1697, ...,  617,  620,  578]],
    
           [[1189,  971,  799, ...,  831,  917,  951],
            [1058,  883,  739, ...,  864,  858,  923],
            [ 894,  764,  684, ...,  895,  819,  919],
            ...,
            [ 972,  942,  903, ...,  695,  565,  538],
            [1207,  927,  865, ...,  820,  803,  700],
            [1983, 1142, 1016, ...,  911,  877,  808]]], dtype=int16)

  • nbart_blue
    (time, y, x)
    int16
    1170 1168 1175 1171 ... 670 685 642

    units :

    1

    nodata :

    -999

    crs :

    EPSG:3577

    grid_mapping :

    spatial_ref

    array([[[1170, 1168, 1175, ..., 1008, 1058, 1020],
            [1171, 1182, 1185, ...,  994, 1029, 1022],
            [1180, 1186, 1205, ...,  955,  968, 1017],
            ...,
            [1472, 1466, 1459, ..., 1066, 1052, 1054],
            [1400, 1407, 1440, ..., 1126, 1077, 1067],
            [1356, 1477, 1482, ..., 1136, 1117, 1111]],
    
           [[3964, 4035, 4169, ..., 5398, 5566, 5292],
            [3629, 3744, 3897, ..., 5256, 5409, 5154],
            [4040, 4120, 4365, ..., 4900, 4901, 4871],
            ...,
            [2790, 2709, 2497, ..., 5580, 5436, 5382],
            [2710, 2627, 2545, ..., 5580, 5430, 5346],
            [2596, 2588, 2555, ..., 5480, 5413, 5434]],
    
           [[ 729,  716,  731, ...,  572,  645,  542],
            [ 699,  712,  722, ...,  554,  576,  541],
            [ 710,  734,  752, ...,  497,  482,  525],
            ...,
    ...
            ...,
            [1155, 1516, 2219, ..., 2827, 3148, 3584],
            [1341, 1791, 2239, ..., 3171, 3285, 3615],
            [1360, 1805, 2264, ..., 3768, 3840, 3770]],
    
           [[ 473,  384,  401, ...,  472,  546,  534],
            [ 419,  399,  389, ...,  445,  461,  466],
            [ 383,  380,  384, ...,  393,  367,  414],
            ...,
            [1214, 1251, 1145, ...,  322,  271,  237],
            [1101, 1001, 1249, ...,  370,  376,  292],
            [1044, 1308, 1286, ...,  405,  408,  403]],
    
           [[ 963,  746,  634, ...,  810,  860,  884],
            [ 842,  687,  594, ...,  828,  801,  884],
            [ 692,  578,  525, ...,  861,  809,  911],
            ...,
            [ 898,  853,  840, ...,  533,  474,  460],
            [1064,  867,  831, ...,  592,  578,  572],
            [1909, 1111,  971, ...,  670,  685,  642]]], dtype=int16)

• Attributes: (2)

  crs :

  EPSG:3577

  grid_mapping :

  spatial_ref

Plot an rgb image to confirm we have data

The white regions are cloud cover.

rgb(ds, index=0, percentile_stretch=(0.1, 0.9))

Export GeoTIFFs

Single-band, single time-slice data

This method uses the datacube.utils.cog function write_cog (where COG stands for Cloud Optimised GeoTIFF) to export a simple single-band, single time-slice GeoTIFF file. A few important caveats should be noted when using this function:

1. It requires an xarray.DataArray; supplying an xarray.Dataset will return an error. To convert an xarray.Dataset to an xarray.DataArray run the following:

da = ds.to_array()

2. This function generates a temporary in-memory GeoTIFF file without compression. This means the function will temporarily use about 1.5 to 2 times the memory of the input xarray.DataArray

# Select a single time-slice and a single band from the dataset.
singleband_da = ds.nbart_red.isel(time=0)

# Write GeoTIFF to a location
write_cog(geo_im=singleband_da,
          fname='red_band.tif',
          overwrite=True)

PosixPath('red_band.tif')

Multi-band, single time-slice data

Here we select a single time and export all the bands in the dataset using write_cog:

# Select a single time-slice
rgb_da = ds.isel(time=0).to_array()

# Write multi-band GeoTIFF to a location
write_cog(geo_im=rgb_da,
          fname='rgb.tif',
          overwrite=True)

PosixPath('rgb.tif')

Multi-band, multiple time-slice data

If we want to export all of the time steps in a dataset, we can wrap write_cog in a for-loop and export each time slice as an individual GeoTIFF file:

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

    # We will use the date of the satellite image to name the GeoTIFF
    date = ds.isel(time=i).time.dt.strftime('%Y-%m-%d').data
    print(f'Writing {date}')

    # Convert current time step into a `xarray.DataArray`
    singletimestamp_da = ds.isel(time=i).to_array()

    # Write GeoTIFF
    write_cog(geo_im=singletimestamp_da,
              fname=f'{date}.tif',
              overwrite=True)

Writing 2020-01-07
Writing 2020-01-16
Writing 2020-01-23
Writing 2020-02-01
Writing 2020-02-08
Writing 2020-02-17
Writing 2020-02-24
Writing 2020-03-04
Writing 2020-03-11
Writing 2020-03-20
Writing 2020-03-27

Exporting GeoTIFFs from a dask array

Note: For more information on using dask, refer to the Parallel processing with Dask notebook

If you pass a lazily-loaded dask array into the function, write_cog will not immediately output a GeoTIFF, but will instead return a dask.delayed object:

# Lazily load data using dask
ds_dask = dc.load(product='ga_ls8c_ard_3',
                  dask_chunks={},
                  **query)

# Run `write_cog`
ds_delayed = write_cog(geo_im=ds_dask.isel(time=0).to_array(),
                       fname='dask_geotiff.tif',
                       overwrite=True)

# Print dask.delayed object
print(ds_delayed)

Delayed('_write_cog-94418718-68c7-4e80-a786-c75e7d6fe50a')

To trigger the GeoTIFF to be exported to file, run .compute() on the dask.delayed object. The file will now appear in the file browser to the left.

ds_delayed.compute()

PosixPath('dask_geotiff.tif')

Advanced

Passing custom rasterio parameters

By default, write_cog will read attributes like nodata directly from the attributes of the data itself. However, write_cog also accepts any parameters that could be passed to the underlying `rasterio.open <https://rasterio.readthedocs.io/en/latest/api/rasterio.io.html>`__ function, so these can be used to override the default attribute values.

For example, it can be useful to provide a new nodata value, for instance if data is transformed to a different dtype or scale and the original nodata value is no longer appropriate. By default, write_cog will use a nodata value of -999 as this is what is stored in the dataset’s attributes:

# Select a single time-slice and a single band from the dataset.
singleband_da = ds.nbart_red.isel(time=0)

# Print nodata attribute value
print(singleband_da.nodata)

-999

To override this value and use a nodata value of 0 instead:

# Write GeoTIFF
write_cog(geo_im=singleband_da,
          fname='custom_nodata.tif',
          overwrite=True,
          nodata=0.0)

PosixPath('custom_nodata.tif')

We can verify the nodata value is now set to 0.0 by reading the data back in with rasterio:

with rasterio.open('custom_nodata.tif') as geotiff:
    print(geotiff.nodata)

0.0

Unsetting nodata for float datatypes containing NaN

A common analysis workflow is to load data from the datacube, then mask out certain values by setting them to NaN. For example, we may use the .where() method to set all -999 nodata values in our data to NaN:

Note: The mask_invalid_data function from datacube.utils.masking can also be used to automatically set nodata values from the data’s attributes to NaN. This function will also automatically drop the nodata attribute from the data, removing the need for the steps below. For more information, see the Masking data notebook.

# Select a single time-slice and a single band from the dataset.
singleband_da = ds.nbart_red.isel(time=0)

# Set -999 values to NaN
singleband_masked = singleband_da.where(singleband_da != -999)
singleband_masked

<xarray.DataArray 'nbart_red' (y: 229, x: 356)>
array([[1691., 1729., 1738., ..., 1180., 1347., 1175.],
       [1680., 1731., 1763., ..., 1140., 1220., 1158.],
       [1737., 1778., 1816., ..., 1070., 1047., 1131.],
       ...,
       [1766., 1808., 1830., ..., 1175., 1069., 1050.],
       [1748., 1731., 1687., ..., 1423., 1193., 1119.],
       [1520., 1802., 1868., ..., 1417., 1340., 1296.]])
Coordinates:
    time         datetime64[ns] 2020-01-07T23:56:40.592891
  * y            (y) float64 -3.955e+06 -3.955e+06 ... -3.962e+06 -3.962e+06
  * x            (x) float64 1.544e+06 1.544e+06 ... 1.554e+06 1.554e+06
    spatial_ref  int32 3577
Attributes:
    units:         1
    nodata:        -999
    crs:           EPSG:3577
    grid_mapping:  spatial_ref

xarray.DataArray

'nbart_red'

• y: 229
• x: 356

• 1.691e+03 1.729e+03 1.738e+03 ... 1.417e+03 1.34e+03 1.296e+03

  array([[1691., 1729., 1738., ..., 1180., 1347., 1175.],
         [1680., 1731., 1763., ..., 1140., 1220., 1158.],
         [1737., 1778., 1816., ..., 1070., 1047., 1131.],
         ...,
         [1766., 1808., 1830., ..., 1175., 1069., 1050.],
         [1748., 1731., 1687., ..., 1423., 1193., 1119.],
         [1520., 1802., 1868., ..., 1417., 1340., 1296.]])

• Coordinates: (4)
  • time
    ()
    datetime64[ns]
    2020-01-07T23:56:40.592891

    units :

    seconds since 1970-01-01 00:00:00

    array('2020-01-07T23:56:40.592891000', dtype='datetime64[ns]')

  • y
    (y)
    float64
    -3.955e+06 ... -3.962e+06

    units :

    metre

    resolution :

    -30.0

    crs :

    EPSG:3577

    array([-3955425., -3955455., -3955485., ..., -3962205., -3962235., -3962265.])

  • x
    (x)
    float64
    1.544e+06 1.544e+06 ... 1.554e+06

    units :

    metre

    resolution :

    30.0

    crs :

    EPSG:3577

    array([1543665., 1543695., 1543725., ..., 1554255., 1554285., 1554315.])

  • spatial_ref
    ()
    int32
    3577

    spatial_ref :

    PROJCS["GDA94 / Australian Albers",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",132],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3577"]]

    grid_mapping_name :

    albers_conical_equal_area

    array(3577, dtype=int32)

• Attributes: (4)

  units :

  1

  nodata :

  -999

  crs :

  EPSG:3577

  grid_mapping :

  spatial_ref

In this case, since we have replaced -999 nodata values with NaN, the data’s nodata attribute is no longer valid. Before we write our data to file, we should therefore remove the nodata attribute from our data:

# Remove nodata attribute
singleband_masked.attrs.pop('nodata')

# Write GeoTIFF
write_cog(geo_im=singleband_masked,
          fname='masked_data.tif',
          overwrite=True)

PosixPath('masked_data.tif')

If we read this GeoTIFF back in with rasterio, we will see that it no longer has a nodata attribute set:

with rasterio.open('masked_data.tif') as geotiff:
    print(geotiff.nodata)

None

---

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 Discord chat 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 2023

Compatible datacube version:

print(datacube.__version__)

1.8.6

Tags

Tags: sandbox compatible, NCI compatible, GeoTIFF, Cloud Optimised GeoTIFF, dc.load, landsat 8, datacube.utils.cog, write_cog, exporting data