Exportation de fichiers GeoTIFF optimisés pour le nuage

Mots-clés données utilisées ; sentinel-2, méthodes de données ; exportation, format des données ; GeoTIFF

Contexte

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).

Un Cloud Optimized GeoTIFF (COG) est un fichier GeoTIFF ordinaire, destiné à être hébergé sur un serveur de fichiers HTTP, avec une organisation interne qui permet des flux de travail plus efficaces sur le cloud.

Description

Ce carnet de notes montre plusieurs façons d’exporter un fichier GeoTIFF :

  1. Exportation d’un GeoTIFF à bande unique et à tranche de temps unique à partir d’un objet xarray chargé par une requête dc.load.

  2. Exportation d’un fichier GeoTIFF multibande et à tranche de temps unique à partir d’un objet xarray chargé par une requête dc.load.

  3. Exportation de plusieurs fichiers GeoTIFF, un pour chaque tranche de temps d’un objet xarray chargé par une requête dc.load.

  4. Exportation de données à partir de tableaux « dask » chargés paresseusement


Pour commencer

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

Chargement des paquets

[1]:
%matplotlib inline

import datacube
from datacube.utils.cog import write_cog

from deafrica_tools.datahandling import load_ard
from deafrica_tools.plotting import rgb

/env/lib/python3.6/site-packages/geopandas/_compat.py:88: UserWarning: The Shapely GEOS version (3.7.2-CAPI-1.11.0 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow.
  shapely_geos_version, geos_capi_version_string

Se connecter au datacube

[2]:
dc = datacube.Datacube(app='Exporting_GeoTIFFs')

Chargement des données Sentinel-2 depuis le datacube

Nous chargeons ici une série chronologique d’images satellite Sentinel-2 via l’API datacube. Cela nous fournira des données avec lesquelles travailler.

[3]:
lat, lon = 13.94, -16.54
buffer = 0.05

# Create a reusable query
query = {
    'x': (lon-buffer, lon+buffer),
    'y': (lat+buffer, lat-buffer),
    'time': ('2020-11-01', '2020-12-15'),
    'resolution': (-20, 20),
    'measurements': ['red', 'green', 'blue'],
    'output_crs':'EPSG:6933'
}

# Load available data from Landsat 8 and filter to retain only times
# with at least 50% good data
ds = load_ard(dc=dc,
              products=['s2_l2a'],
              **query)

# Print output data
print(ds)

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Applying pixel quality/cloud mask
Loading 9 time steps
<xarray.Dataset>
Dimensions:      (time: 9, x: 483, y: 620)
Coordinates:
  * time         (time) datetime64[ns] 2020-11-03T11:47:42 ... 2020-12-13T11:...
  * y            (y) float64 1.768e+06 1.768e+06 ... 1.755e+06 1.755e+06
  * x            (x) float64 -1.601e+06 -1.601e+06 ... -1.591e+06 -1.591e+06
    spatial_ref  int32 6933
Data variables:
    red          (time, y, x) float32 nan nan nan nan ... 110.0 131.0 140.0
    green        (time, y, x) float32 nan nan nan nan ... 247.0 248.0 290.0
    blue         (time, y, x) float32 nan nan nan nan ... 99.0 67.0 87.0 89.0
Attributes:
    crs:           EPSG:6933
    grid_mapping:  spatial_ref

Tracez une image rgb pour confirmer que nous avons des données.

Les régions blanches correspondent à la couverture nuageuse.

[4]:
rgb(ds, index=3)

../../../_images/sandbox_notebooks_Frequently_used_code_Exporting_GeoTIFFs_12_0.png

Exportation d’une bande unique, d’une tranche de temps unique GeoTIFF

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 COG.

Quelques mises en garde importantes sont à noter lors de l’utilisation de cette fonction : 1. Elle nécessite un xarray.DataArray ; fournir un xarray.Dataset retournera une erreur. Pour convertir un xarray.Dataset en un tableau, exécutez ce qui suit :

da = ds.to_array()
  1. Cette fonction génère un fichier tiff temporaire en mémoire, sans compression. Cela signifie que la fonction utilisera environ 1,5 à 2 fois la mémoire nécessaire en utilisant la fonction dépréciée datacube.helper.write_geotiff.

  2. Si vous passez un dask array dans la fonction, write_cog ne produira pas de geotiff, mais retournera un objet Dask Delayed. Pour déclencher la sortie du geotiff, exécutez .compute() sur l’objet Dask Delayed :

    write_cog(ds.red.isel(time=0), "red.tif").compute()
    
[5]:
# Select a single time-slice and a single band from the dataset.
singleBandtiff = ds.red.isel(time=5)

# Write GeoTIFF to a location
write_cog(singleBandtiff,
          fname="red_band.tif",
          overwrite=True)

[5]:
PosixPath('red_band.tif')

Exportation d’un GeoTIFF multi-bandes, à tranche de temps unique.

Ici, nous sélectionnons un seul temps et exportons toutes les bandes de l’ensemble de données en utilisant la fonction datacube.utils.cog.write_cog.

[6]:
# Select a single time-slice
rgb_tiff = ds.isel(time=5).to_array()

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

[6]:
PosixPath('rgb.tif')

Exportation de plusieurs GeoTIFFs, un pour chaque tranche de temps d’un xarray.

Si nous voulons exporter tous les pas de temps d’un ensemble de données comme un GeoTIFF, nous pouvons envelopper notre fonction write_cog dans une boucle for.

[7]:
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 = ds.isel(time=i).to_array()

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

Writing 2020-11-03
Writing 2020-11-08
Writing 2020-11-13
Writing 2020-11-18
Writing 2020-11-23
Writing 2020-11-28
Writing 2020-12-03
Writing 2020-12-08
Writing 2020-12-13

Avancé

Exportation de GeoTIFFs à partir d’un tableau dask

Si vous passez un tableau dask chargé paresseusement dans la fonction, write_cog ne produira pas immédiatement un GeoTIFF, mais retournera un objet dask.delayed :

[8]:
# Lazily load data using dask
ds_dask = dc.load(product='s2_l2a',
                  dask_chunks={},
                  **query)

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

# Print dask.delayed object
print(ds_delayed)
Delayed('_write_cog-54a5f520-6745-4265-bd08-10031f66d67f')

Pour déclencher l’exportation du GeoTIFF vers un fichier, exécutez .compute() sur l’objet dask.delayed. Le fichier apparaîtra maintenant dans le navigateur de fichiers à gauche.

[9]:
ds_delayed.compute()
[9]:
PosixPath('dask_geotiff.tif')

Informations complémentaires

Licence: Le code de ce carnet est sous licence Apache License, Version 2.0. Les données de Digital Earth Africa sont sous licence Creative Commons by Attribution 4.0.

Contact: Si vous avez besoin d’aide, veuillez poster une question sur le canal Slack Open Data Cube ou sur le GIS Stack Exchange en utilisant la balise open-data-cube (vous pouvez voir les questions précédemment posées ici). Si vous souhaitez signaler un problème concernant ce carnet de notes, vous pouvez en déposer un sur Github.

Version du datacube compatible:

[10]:
print(datacube.__version__)
1.8.4.dev52+g07bc51a5

Dernier test:

[11]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[11]:
'2021-04-15'