Reprojection de données de cube de données et de données raster

• Produits utilisés : gm_s2_annual

• Exigences spéciales : ce bloc-notes charge les données d’un fichier raster externe (madagascar_CHPclim_02.tif) à partir du dossier Supplementary_data de ce référentiel

Mots clés données utilisées ; géomédiane sentinel-2, données utilisées ; CHIRPS, méthodes de données ; reprojection, méthodes de données ; rééchantillonnage, format de données ; GeoTIFF, méthodes de données ; geobox

Aperçu

Il est souvent utile de combiner les données du datacube avec d’autres jeux de données raster externes. Par exemple, nous pouvons vouloir utiliser un raster de précipitations pour concentrer notre analyse sur les données satellite afin de masquer les zones de précipitations plus ou moins élevées. Cependant, il peut être difficile de combiner des jeux de données s’ils ont des étendues, des résolutions (par exemple 20 m contre 250 m) ou des systèmes de référence de coordonnées différents (par exemple WGS 84 contre Africa Albers Equal Area Conic). Pour pouvoir combiner ces jeux de données, nous devons pouvoir les reprojeter dans des grilles spatiales identiques avant l’analyse.

Datacube stocke des informations sur l’emplacement, la résolution et le système de référence de coordonnées d’une grille rectangulaire de données à l’aide d’un objet appelé GeoBox. Cet objet GeoBox` est ajouté à toutes les données chargées à partir du datacube et à tous les rasters chargés avec xr.open_rasterio (à condition que import datacube soit exécuté avant que le raster ne soit chargé). Les fonctions Datacube peuvent utiliser cet objet pour fournir un modèle qui peut être utilisé pour reprojeter les données raster et datacube - soit en appliquant cette reprojection directement lorsque de nouvelles données sont chargées, soit pour reprojeter les données existantes qui ont déjà été chargées en mémoire. Cela permet de reprojeter facilement un ensemble de données dans la grille spatiale exacte d’un autre ensemble de données pour une analyse plus approfondie.

Description

Ce bloc-notes montre comment exécuter trois flux de travail de reprojection clés lorsque vous travaillez avec des données de cube de données et des rasters externes :

1. Chargement et reprojection des données du Datacube pour qu’elles correspondent à un fichier raster à l’aide de « dc.load »

2. Reprojection des données de datacube existantes pour correspondre à un raster à l’aide de « xr_reproject »

3. Chargement et reprojection d’un raster pour correspondre aux données du Datacube à l’aide de rio_slurp_xarray

---

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

Importez les packages Python utilisés pour l’analyse.

import datacube
import rioxarray
import matplotlib.pyplot as plt
from odc.algo import xr_reproject
from datacube.testutils.io import rio_slurp_xarray
from datacube.utils.masking import mask_invalid_data

from deafrica_tools.plotting import rgb

Se connecter au datacube

Connectez-vous au datacube pour que nous puissions accéder aux données de DE Africa. Le paramètre « app » est un nom unique pour l’analyse qui est basé sur le nom du fichier du notebook.

dc = datacube.Datacube(app="Reprojecting_data")

Chargement et reprojection des données du Datacube pour correspondre à un raster

Charger un fichier raster

Nous chargeons d’abord un raster GeoTIFF à partir d’un fichier en utilisant « rioxarray.open_rasterio ». L’exemple ci-dessous charge des données de précipitations moyennes à long terme pour le mois de février conformément à « Climate Hazards Group InfraRed Precipitation with Station <https://www.chc.ucsb.edu/data/chirps> » (CHIRPS) pour une section de Madagascar. Il a une résolution spatiale d’environ 0,05 degré décimal (environ 5 km de pixels à l’équateur) dans le système de référence de coordonnées WGS 84 (EPSG:4326).

# Path to raster file
raster_path = "../Supplementary_data/Reprojecting_data/madagascar_CHPclim_02.tif"

# Load raster, and remove redundant "band" dimension
raster = rioxarray.open_rasterio(raster_path).squeeze("band", drop=True)

print(raster)

<xarray.DataArray (y: 20, x: 20)>
[400 values with dtype=float32]
Coordinates:
  * x            (x) float64 46.03 46.08 46.13 46.18 ... 46.83 46.88 46.93 46.98
  * y            (y) float64 -18.03 -18.08 -18.13 ... -18.88 -18.93 -18.98
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

Si nous traçons notre raster chargé, nous pouvons voir qu’il y a des zones avec des précipitations plus élevées et des zones avec des précipitations plus faibles.

raster.plot(size=6)

<matplotlib.collections.QuadMesh at 0x7fa4bdd79e70>

Objets GeoBox

Now we have loaded our raster dataset, we can inspect its GeoBox object that we will use to allow us to reproject data. The GeoBox can be accessed using the .geobox method. It includes a set of information that together completely define the spatial grid of our data:

• The width (e.g. 20) and height (e.g. 20) of our data in pixels

• An Affine object which defines the spatial resolution (e.g. -0.05000000074505806 and 0.05000000074505806) and spatial position (e.g. 46.00000336766243 and -18.00000160932541) of our data

• The coordinate reference system of our data (e.g. +init=epsg:4326)

raster.geobox

GeoBox(20, 20, Affine(0.05000000074505806, 0.0, 46.00000336766243,
       0.0, -0.05000000074505806, -18.00000160932541), GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]])

La « GeoBox » dispose également d’un certain nombre de méthodes utiles qui peuvent être utilisées pour visualiser des informations sur la grille spatiale de nos données. Par exemple, nous pouvons inspecter la résolution spatiale de nos données :

raster.geobox.resolution

(-0.05000000074505806, 0.05000000074505806)

Ou obtenir des informations sur l’étendue spatiale des données :

raster.geobox.extent.boundingbox

BoundingBox(left=46.00000336766243, bottom=-19.00000162422657, right=47.00000338256359, top=-18.00000160932541)

Remarque : pour plus d’informations sur les objets « GeoBox » et une liste complète de leurs méthodes, reportez-vous à la documentation de datacube <https://datacube-core.readthedocs.io/en/latest/api/core-classes/generate/datacube.utils.geometry.GeoBox.html>`__.

Charger et reprojeter les données du Datacube

Nous pouvons maintenant utiliser Datacube pour charger et reprojeter des données satellite afin qu’elles correspondent exactement au système de référence de coordonnées et à la résolution de nos données raster. En spécifiant « like=raster.geobox », nous pouvons charger des données Datacube qui seront reprojetées pour correspondre à la grille spatiale de notre raster.

# Load and reproject data from datacube
ds = dc.load(product="gm_s2_annual",
             measurements=["red", "green", "blue"],
             time=("2018"),
             like=raster.geobox,
             resampling="nearest",
             group_by="solar_day")

Lorsque nous traçons le résultat, il doit apparaître de manière similaire au raster que nous avons chargé ci-dessus, qui a une résolution inférieure.

rgb(ds, col='time')

Nous pouvons également comparer directement la géobox des deux ensembles de données pour vérifier qu’ils partagent la même grille spatiale :

ds.geobox == raster.geobox

True

Maintenant que nos deux jeux de données partagent la même grille spatiale, nous pouvons utiliser notre raster comme masque. Par exemple, nous pouvons masquer tous les pixels satellites à l’exception de ceux présentant des précipitations plus élevées (par exemple, plus de 350 mm/mois de pluie) :

# Rename raster dimensions to match datacube conventions for data
# with geographic coordinates
if raster.geobox.crs.geographic:
    raster = raster.rename({"x": "longitude", "y": "latitude"})

# Identify high rainfall areas
high_rainfall = raster > 350

# Apply mask to set lower-rainfall areas to `NaN`
ds_masked = ds.where(high_rainfall)

# Plot the masked data
rgb(ds_masked)

/usr/local/lib/python3.10/dist-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

Reprojection des données de datacube existantes pour correspondre à un raster

L’exemple ci-dessus montre comment charger de nouvelles données satellite à partir du datacube pour correspondre à la grille spatiale d’un raster. Cependant, il se peut que nous ayons déjà chargé des données satellite avec un système de référence de coordonnées et une résolution différents de ceux de notre raster. Dans ce cas, nous devons reprojeter cet ensemble de données existant pour qu’il corresponde à notre raster.

Par exemple, nous avons chargé les données géomédianes Sentinel-2 2018 à partir du datacube avec des pixels de résolution 20 m dans le système de référence de coordonnées « EPSG:6933 ». Notez que dans cet exemple, nous spécifions manuellement les paramètres « x », « y », « resolution » et « output_crs », plutôt que de les prendre directement à partir de notre raster en utilisant « like=raster.geobox » dans l’exemple précédent.

# Load data from datacube
ds = dc.load(product="gm_s2_annual",
             measurements=["red", "green", "blue"],
             time=("2018"),
             x=(46.00, 47.00),
             y=(-19, -18),
             resolution=(-20, 20),
             output_crs="EPSG:6933",
             group_by="solar_day")

# Load raster, and remove redundant "band" dimension
raster = rioxarray.open_rasterio(raster_path).squeeze("band", drop=True)

Si nous traçons nos données satellite, nous pouvons voir qu’elles ont une résolution bien supérieure à celle de notre raster de précipitations CHIRPS pixellisé :

rgb(ds)

Reprojeter les données du datacube

Nous pouvons maintenant utiliser la fonction « xr_reproject » pour reprojeter notre jeu de données satellite haute résolution existant. Nous spécifions « geobox=raster.geobox » pour demander que les données soient reprojetées pour correspondre à la grille spatiale de notre raster basse résolution.

Pour contrôler la manière dont les données sont reprojetées, nous pouvons spécifier une méthode de « rééchantillonnage » personnalisée qui contrôlera la manière dont nos pixels haute résolution seront transformés en pixels de résolution inférieure. Dans ce cas, nous spécifions « moyenne », qui définira la valeur de chaque pixel plus grand sur la moyenne de tous les pixels plus petits qui se trouvent dans sa limite de pixels.

Remarque : reportez-vous à la « documentation rasterio <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio-enums-module> » pour obtenir la liste complète des méthodes de rééchantillonnage disponibles.

# Temporary workaround for bug in `xr_reproject`
if raster.geobox.crs.geographic:
    ds = ds.rename({"x": "longitude", "y": "latitude"})

# Reproject data
ds_reprojected = xr_reproject(src=ds,
                              geobox=raster.geobox,
                              resampling="average")

#Set nodata to `NaN`
ds_reprojected = mask_invalid_data(ds_reprojected)

Maintenant, si nous traçons notre ensemble de données reprojetées, nous pouvons voir que notre imagerie satellite a maintenant une résolution similaire à notre raster basse résolution (par exemple avec une apparence pixellisée). Notez cependant que ce résultat apparaîtra plus lisse que l’exemple précédent en raison de la méthode de rééchantillonnage « moyenne » utilisée ici (par rapport à « la plus proche »).

rgb(ds_reprojected)

Une fois de plus, nous pouvons également vérifier que les deux ensembles de données ont des grilles spatiales identiques :

ds_reprojected.geobox == raster.geobox

True

Chargement et reprojection d’un raster pour correspondre aux données du Datacube

Au lieu de reprojeter les données satellite pour qu’elles correspondent à la résolution et au système de projection de notre raster, nous pouvons souhaiter reprojeter notre raster pour qu’il corresponde à la grille spatiale de nos données satellite. Cela peut être particulièrement utile lorsque nous avons un fichier raster à faible résolution (par exemple, comme les données CHIRPS d’une résolution d’environ 5 km que nous utilisons dans cet exemple), mais que nous ne voulons pas perdre la meilleure résolution spatiale de nos données satellite.

Charger les données du Datacube

Comme dans l’exemple précédent, nous pouvons charger des données satellite à partir du datacube avec une résolution spatiale de 20 m et un système de référence de coordonnées « EPSG:6933 » :

# Load data from datacube
ds = dc.load(product="gm_s2_annual",
             measurements=["red", "green", "blue"],
             time=("2018"),
             x=(46.00, 47.00),
             y=(-19, -18),
             resolution=(-20, 20),
             output_crs="EPSG:6933",
             group_by="solar_day")

# Plot 20 m resolution satellite data
rgb(ds)

Charger et reprojeter des données raster

Nous pouvons maintenant utiliser la fonction rio_slurp_xarray pour charger et reprojeter notre fichier raster afin qu’il corresponde à notre jeu de données satellite à plus haute résolution. Nous spécifions gbox=ds.geobox pour demander que notre raster soit reprojeté pour correspondre à la grille spatiale de ds. Comme dans l’exemple précédent de xr_reproject, nous pouvons également spécifier une méthode de resampling personnalisée qui sera utilisée pendant la reprojection. Dans ce cas, nous spécifions 'bilinéaire', qui produira une sortie lisse sans limites de pixels évidentes.

Remarque : reportez-vous à la « documentation rasterio <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio-enums-module> » pour obtenir la liste complète des méthodes de rééchantillonnage disponibles.

# Load raster and reproject to match satellite dataset
raster_reprojected = rio_slurp_xarray(fname=raster_path,
                                      gbox=ds.geobox,
                                      resampling="bilinear")

# Set nodata to `NaN`
raster_reprojected = mask_invalid_data(raster_reprojected)

Si nous traçons nos données raster rééchantillonnées, elles devraient maintenant apparaître moins pixelisées que le raster d’origine d’une résolution d’environ 250 m.

raster_reprojected.plot(size=6, vmin=200);

Le raster rééchantillonné doit également correspondre à la grille spatiale de nos données satellite à plus haute résolution :

raster_reprojected.geobox == ds.geobox

True

Maintenant que nos deux jeux de données partagent la même grille spatiale, nous pouvons utiliser notre raster rééchantillonné pour masquer notre jeu de données satellite à plus haute résolution comme nous l’avons fait dans la première section (par exemple, masquer toutes les zones avec moins de 350 mm de précipitations en février). Par rapport au premier exemple, ce jeu de données satellite masqué devrait apparaître avec une résolution beaucoup plus élevée, avec une pixellisation beaucoup moins évidente :

# Identify high rainfall areas
high_rainfall = raster_reprojected > 350

# Apply mask to set high rainfall areas to `NaN`
ds_masked = ds.where(high_rainfall)

# Plot the masked data
rgb(ds_masked)

/usr/local/lib/python3.10/dist-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

---

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 :

print(datacube.__version__)

1.8.15

Dernier test :

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

'2023-08-14'