Exportation de fichiers GeoTIFF optimisés pour le cloud

  • Produits utilisés : s2_l2a

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

Aperçu

À la fin d’une analyse, il peut être utile d’exporter les données vers un fichier GeoTIFF (par exemple outputname.tif), soit pour enregistrer les résultats, soit pour permettre l’exploration des résultats dans une plateforme logicielle SIG (par exemple ArcGIS ou QGIS).

Un « Cloud Optimized GeoTIFF (COG) » est un fichier GeoTIFF standard, 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 bloc-notes présente plusieurs façons d’exporter un fichier GeoTIFF :

  1. Exportation d’un GeoTIFF monobande et monotranche de temps à partir d’un objet xarray chargé via une requête « dc.load »

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

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

  4. Exportation de données à partir de « tableaux dask » chargés de manière différée

Commencer

Pour exécuter cette analyse, exécutez toutes les cellules du bloc-notes, en commençant par la cellule « Charger les packages ».

Charger des paquets

[1]:

%matplotlib inline

import datacube
import geopandas as gpd
from datacube.utils.cog import write_cog
from datacube.utils.geometry import Geometry

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

Se connecter au datacube

[2]:

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

Charger les données Sentinel-2 à partir du 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.

La cellule suivante définit les paramètres qui définissent la zone d’intérêt et la durée de l’analyse. Les paramètres sont

  • « lat » : la latitude centrale à analyser (par exemple « 6,502 »).

  • « lon » : la longitude centrale à analyser (par exemple « -1,409 »).

  • « buffer » : le nombre de degrés carrés à charger autour de la latitude et de la longitude centrales. Pour des temps de chargement raisonnables, définissez cette valeur sur « 0,1 » ou moins.

Pour définir la zone d’intérêt, deux méthodes sont disponibles :

  1. En spécifiant la latitude, la longitude et la zone tampon. Cette méthode nécessite que vous saisissiez la latitude centrale, la longitude centrale et la valeur de la zone tampon en degrés carrés autour du point central que vous souhaitez analyser. Par exemple, « lat = 10,338 », « lon = -1,055 » et « buffer = 0,1 » sélectionneront une zone avec un rayon de 0,1 degré carré autour du point avec les coordonnées (10,338, -1,055).

  2. By uploading a polygon as a GeoJSON or Esri Shapefile. If you choose this option, you will need to upload the geojson or ESRI shapefile into the Sandbox using Upload Files button c78e43a9b3594be7870e6d98f9b77bf3 in the top left corner of the Jupyter Notebook interface. ESRI shapefiles must be uploaded with all the related files (.cpg, .dbf, .shp, .shx). Once uploaded, you can use the shapefile or geojson to define the area of interest. Remember to update the code to call the file you have uploaded.

Pour utiliser l’une de ces méthodes, vous pouvez décommenter la ligne de code concernée et commenter l’autre. Pour commenter une ligne, ajoutez le symbole "#" avant le code que vous souhaitez commenter. Par défaut, la première option qui définit l’emplacement à l’aide de la latitude, de la longitude et du tampon est utilisée.

[3]:

# Set the area of interest

# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=13.94, lon=-16.54, buffer=0.05)

# Method 2: Use a polygon as a GeoJSON or Esri Shapefile.
# aoi = define_area(vector_path='aoi.shp')

#Create a geopolygon and geodataframe of the area of interest
geopolygon = Geometry(aoi["features"][0]["geometry"], crs="epsg:4326")
geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)

# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])

# Create a reusable query
query = {
    'x': lon_range,
    'y': lat_range,
    '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

/usr/local/lib/python3.10/dist-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(

<xarray.Dataset>
Dimensions:      (time: 9, y: 620, x: 483)
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 RVB pour confirmer que nous avons des données

Les régions blanches sont une couverture nuageuse.

[4]:

rgb(ds, index=3)
../../../_images/sandbox_notebooks_Frequently_used_code_Exporting_GeoTIFFs_12_0.png

Exporter un GeoTIFF à bande unique et à tranche de temps unique

Cette méthode utilise la fonction datacube.utils.cog write_cog, où cog signifie Cloud-Optimised-Geotiff pour exporter un COG simple à bande unique et à tranche de temps unique.

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 a xarray.Dataset to an array run the following:

    da = ds.to_array()

  2. 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 requise en utilisant la fonction dépréciée « datacube.helper.write_geotiff ».

  3. Si vous passez un « tableau dask » dans la fonction, « write_cog » ne générera pas de géotiff, mais renverra à la place un objet « Dask Delayed ». Pour déclencher la sortie du géotiff, exécutez « .compute() » sur l’objet dask delay :

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

Exporter un GeoTIFF multibande à tranche temporelle unique

Ici, nous sélectionnons une seule heure et exportons toutes les bandes de l’ensemble de données à l’aide de 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')

Exporter plusieurs GeoTIFF, un pour chaque tranche de temps d’un xarray

Si nous voulons exporter tous les pas de temps d’un ensemble de données au format 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 fichiers GeoTIFF à partir d’une matrice Dask

Si vous passez un tableau dask chargé paresseusement dans la fonction, write_cog ne générera pas immédiatement un GeoTIFF, mais renverra à la place 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-5ff4ca46-32b4-44c9-9529-5cfea7c9b1a3')

Pour déclencher l’exportation du GeoTIFF vers un fichier, exécutez .compute() sur l’objet dask.delayed. Le fichier apparaîtra alors 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, version 2.0 <https://www.apache.org/licenses/LICENSE-2.0>. Les données de Digital Earth Africa sont sous licence Creative Commons par attribution 4.0 <https://creativecommons.org/licenses/by/4.0/>.

Contact : Si vous avez besoin d’aide, veuillez poster une question sur le canal Slack Open Data Cube <http://slack.opendatacube.org/>`__ ou sur le GIS Stack Exchange en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment ici). Si vous souhaitez signaler un problème avec ce bloc-notes, vous pouvez en déposer un sur Github.

Version de Datacube compatible :

[10]:

print(datacube.__version__)

1.8.15

Dernier test :

[11]:

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

[11]:

'2023-08-14'