Travailler avec le temps dans xarray

Produits utilisés : s2_l2a

Mots clés analyse; séries chronologiques, données utilisées; sentinel-2, méthodes de données; groupby,:index:méthodes de données; le plus proche, index:méthodes de données; interpolation, méthodes de données; rééchantillonnage, méthodes de données; composition

Aperçu

Les données de séries chronologiques <https://en.wikipedia.org/wiki/Time_series> sont une série de points de données généralement capturés à des moments successifs dans le temps. Dans un contexte de télédétection, les données de séries chronologiques sont une séquence d’images satellite discrètes prises dans la même zone à des moments successifs. L’analyse des séries chronologiques utilise différentes méthodes pour extraire des statistiques, des modèles et d’autres caractéristiques significatives des données. Les données et l’analyse des séries chronologiques ont de nombreuses applications allant de la surveillance des cultures agricoles à la détection des changements naturels de la végétation, en passant par la cartographie de la prospectivité minérale et la modélisation de la hauteur des marées.

Description

Le package Python « xarray » fournit de nombreuses techniques utiles pour traiter les données de séries chronologiques qui peuvent être appliquées aux données de Digital Earth Africa. Ce bloc-notes montre comment utiliser les techniques « xarray » pour :

Sélectionnez différentes périodes de données (par exemple, année, mois, jour) à partir d’un « xarray.Dataset »
Utiliser les accesseurs datetime pour extraire des informations supplémentaires de la dimension « temps » d’un ensemble de données
Résumer les données de séries chronologiques pour différentes périodes à l’aide de « .groupby() » et « .resample() »
Interpoler des données de séries chronologiques pour estimer les conditions du paysage à une date précise que le satellite n’a pas observées

Pour plus d’informations sur les techniques présentées ci-dessous, reportez-vous au guide des données de séries chronologiques « xarray » <http://xarray.pydata.org/en/stable/time-series.html>`__.

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 matplotlib.pyplot as plt
import numpy as np
import geopandas as gpd
from datacube.utils.geometry import Geometry

from deafrica_tools.datahandling import load_ard, mostcommon_crs
from deafrica_tools.areaofinterest import define_area

Se connecter au datacube

[2]:

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

Chargement des données Landsat

Tout d’abord, nous chargeons environ deux ans de données Sentinel-2, en utilisant la fonction load_ard et en filtrant les pas de temps avec au moins 95 % de pixels de bonne qualité.

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

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

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

# 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': ('2018-01', '2019-12'),
    'resolution': (-20, 20),
    'measurements':['red', 'green', 'blue', 'nir']
}

# Identify the most common projection system in the input query
output_crs = mostcommon_crs(dc=dc, product='s2_l2a', query=query)

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

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Counting good quality pixels for each time step

/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(
/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(

Filtering to 34 out of 144 time steps with at least 95.0% good quality pixels
Applying pixel quality/cloud mask
Loading 34 time steps

Explorer les données de matrice X en utilisant le temps

Nous allons ici explorer plusieurs façons d’utiliser la dimension temporelle dans un « xarray.Dataset ». Cette section décrit la sélection, la synthèse et l’interpolation des données à des moments précis.

Regroupement et rééchantillonnage par temps

« xarray » fournit également des raccourcis pour agréger des données au fil du temps. Dans l’exemple ci-dessous, nous regroupons d’abord nos données par saison, puis prenons la médiane de chaque groupe. Cela produit un nouvel ensemble de données avec seulement quatre observations (une par saison).

[13]:

# Group the time series into seasons, and take median of each time period
ds_seasonal = ds.groupby('time.season').median(dim='time')

# Plot the output
ds_seasonal.nir.plot(col='season', col_wrap=4)
plt.show()

../../../_images/sandbox_notebooks_Frequently_used_code_Working_with_time_31_0.png

Nous pouvons également utiliser la méthode .resample() pour résumer notre ensemble de données en blocs de temps plus grands. Dans l’exemple ci-dessous, nous produisons une médiane composite pour chaque 6 mois de données dans notre ensemble de données :

[14]:

# Resample to combine each 6 months of data into a median composite
ds_resampled = ds.resample(time="6m").median()

# Plot the new resampled data
ds_resampled.nir.plot(col="time")
plt.show()

../../../_images/sandbox_notebooks_Frequently_used_code_Working_with_time_33_0.png

Interpolation de nouveaux pas de temps

Parfois, nous souhaitons obtenir des données pour des heures/dates spécifiques qui n’ont pas été observées par un satellite. Pour estimer à quoi ressemblait le paysage à certaines dates, nous pouvons utiliser la méthode .interp() pour interpoler entre les deux observations les plus proches.

Par défaut, la méthode interp() utilise une interpolation linéaire (method='linear'). Une autre option utile est method='nearest', qui renverra l’observation satellite la plus proche de la ou des dates spécifiées.

[15]:

# New dates to interpolate data for
new_dates = ['2018-07-25', '2018-09-01', '2018-12-05']

# Interpolate Landsat values for three new dates
ds_interp = ds.interp(time=new_dates)

# Plot the new interpolated data
ds_interp.nir.plot(col='time')
plt.show()

../../../_images/sandbox_notebooks_Frequently_used_code_Working_with_time_35_0.png

Informations Complémentaires

Licence:

Le code de ce bloc-notes est sous licence « Apache License, Version 2.0 <https://www.apache.org/licenses/LICENSE-2.0> ».

Les données de Digital Earth Africa sont sous licence Creative Commons 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 :

[16]:

print(datacube.__version__)

1.8.15

Dernier test :

[17]:

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

[17]:

'2023-08-14'