In English

ECCC logo

TdM > Aperçu de l'utilisation > Requêtes OGC API - Features avec Python (version script)

Cas d'utilisation : Récupération et affichage de données hydrométriques

Introduction

Les données hydrométriques du Service météorologique du Canada (SMC) d'Environnement et Changement climatique Canada (ECCC) peuvent être facilement accédées via GeoMet-OGC-API. Les services web OGC API - Features disponibles dans GeoMet-OGC-API permettent d'effectuer des requêtes de données géospatiales vectorielles et de récupérer ces données en format GeoJSON. Les fichiers GeoJSON peuvent ensuite être utilisés directement pour réaliser des analyses et afficher les données sur des cartes ou peuvent être soumis à des étapes supplémentaires de traitement de données. Ce cas d'utilisation démontre comment:

  • Effectuer des requêtes et récupérer des données géospatiales vectorielles disponibles dans GeoMet-OGC-API à l'aide de Python
  • Afficher des données de séries temporelles en format graphique et tabulaire
  • Créer une carte interactive présentant les données géospatiales récupérées

Ce cas d'utilisation est disponible en deux versions. Cette version-ci a été concue pour être utilisé de manière traditionnelle dans un terminal ou un environnement de développement intégré, mais peut également être utilisé dans Jupyter Notebook. Si vous recherchez plus d'éléments interactifs tels que des menus déroulants pour sélectionner la valeur des variables, un graphique interactif comportant une plus grande variété de fonctions et une carte interactive basée sur Leaflet, la version interactive de ce cas d'utilisation conçue spécialement pour Jupyter Notebook pourrait vous intéresser.

Une version exécutable de ce cas d'utilisation est disponible. Dans Jupyter Notebook, cliquez sur le symbole ⏩️ ou sélectionnez l'option Cell -> Run All dans la barre de menu pour exécuter l'entièreté de ce carnet. Si les annotations ne sont pas visibles lors du passage de la souris sur le graphique ou sur la carte, exécutez à nouveau la cellule contenant le code du graphique ou de la carte en cliquant dans la cellule puis en cliquant sur le bouton Exécuter dans la barre de menu.

Note importante: Hors de Jupyter Notebook, la commande magique %matplotlib notebook située à la fin de la section du code portant sur les paramètres doit être retirée.

badge

Création d'un outil pour suivre l'évolution des niveaux d'eau

Cathy dirige une compagnie d'excursions en bateau et de sports nautiques à Chilliwack en Colombie-Britannique. Elle aimerait suivre le niveau d'eau de certaines stations hydrométriques situés près de Chilliwack afin de guider la gestion des activités de sa compagnie. Elle aimerait obtenir:

  • Les moyennes journalières de niveau d'eau pour une période de trois mois pour toutes les stations hydrométriques situées dans un rayon de 100 km de Chilliwack
  • Une vue graphique et tabulaire des moyennes journalières de niveau d'eau de cette période
  • Une carte interactive présentant les stations hydrométriques situées dans un rayon de 100 km de Chilliwack et la moyenne journalière de niveau d'eau la plus récente de ces stations pour la période temporelle sélectionnée

Pour y arriver, la première étape est d'importer les modules Python nécessaires et de sélectionner les paramètres de requête désirés.

In [1]:
# Importation des modules
from datetime import date
import json
import math
from textwrap import fill

import cartopy.crs as ccrs
from cartopy.io.img_tiles import OSM
from matplotlib import pyplot as plt, dates as mdates
from osgeo import ogr, osr
from owslib.ogcapi.features import Features
import numpy as np
import pandas as pd
from tabulate import tabulate
In [2]:
# Paramètres

# Coordonnées de Chilliwack
lat = 49.162676
long = -121.958943

# Taille de la zone tampon en kilomètres
buffer = 100

# Dates de début et de fin de la période temporelle pour
# laquelle des données seront récupérées
start_date = date(2018, 6, 1)
end_date = date(2018, 8, 31)

# Code ESPG de la projection désirée pour crééer la zone tampon
# NAD83 / Statistiques Canada Lambert
projection = 3347

# Commentez la ligne suivante si vous utilisez le code hors de Jupyter Notebook
%matplotlib notebook
In [3]:
# Formattage des paramètres pour la requête OGC API - Features

# Zone de délimitation légèrement plus grande que la taille de la zone tampon

# La taille de la zone tampon doit être transformée en degrés afin d'obtenir
# les coordonnées des coins de la zone de délimitation:
# Latitude: 1 km ≈ 0.009° 
# Longitude (au 49e parallèle): 1 km ≈ 0.014°
bbox = [
    str(long - buffer * 0.02),
    str(lat - buffer * 0.01),
    str(long + buffer * 0.02),
    str(lat + buffer * 0.01),
]

# Formattage de la période temporelle sélectionnée
time = f"{start_date}/{end_date}"

Ensuite, les données des stations hydrométriques disponibles sur GeoMet-OGC-API peuvent être récupérées à l'aide d'OWSLib. Comme l'utilisation d'une zone de délimitation à l'aide du paramètre bbox n'est pas aussi précise que l'utilisation d'une zone tampon ronde, GDAL peut être utilisé pour sélectionner uniquement les stations hydrométriques situées strictement dans un rayon de 100 km de Chilliwack.

In [4]:
# Récupération des données de stations hydrométriques
oafeat = Features("https://api.weather.gc.ca/")
station_data = oafeat.collection_items(
    "hydrometric-stations", bbox=bbox, STATUS_EN="Active"
)

# Vérification des données récupérées
if "features" in station_data.keys() and station_data["features"]:
    station_data = json.dumps(station_data, indent=4)
else:
    raise ValueError(
        "Aucune station hydrométrique n'a été trouvée. "
        + "Veuillez revérifier les coordonnées."
    )
In [5]:
# Liste des stations situées dans la zone tampon

# Accéder à la couche de stations hydrométriques
driver = ogr.GetDriverByName("GeoJSON")
data_source = driver.Open(station_data, 0)
layer = data_source.GetLayer()

# Identification du système de référence spatiale (SRS) d'entrée
SRS_input = layer.GetSpatialRef()
SR = osr.SpatialReference(str(SRS_input))
epsg = SR.GetAuthorityCode(None)
SRS_input.ImportFromEPSG(int(epsg))

# Définition du SRS utilisé pour projeter les données
SRS_projected = osr.SpatialReference()
SRS_projected.ImportFromEPSG(projection)

# Transformation du SRS d'entrée à la projection désirée
transform = osr.CoordinateTransformation(SRS_input, SRS_projected)

# Création de la zone tampon pour sélectionner les stations
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(long, lat)
point.Transform(transform)
point_buffer = point.Buffer(buffer * 1000)  # La valeur doit être en mètres

# Sélection des stations situées dans la zone tampon
stations = []

for feature in layer:
    geom = feature.GetGeometryRef().Clone()
    geom.Transform(transform)
    if geom.Intersects(point_buffer):
        stations.append(feature.STATION_NUMBER)

# Générer une erreur si aucune station n'est trouvée
if not stations:
    raise ValueError(
        f"Aucune station hydrométrique n'est située dans un rayon de {buffer}"
        + " km des coordonnées choisies. Veuillez vérifier les coordonnées."
    )

Une fois les stations hydrométriques situées dans un rayon de 100 km de Chilliwack identifiées, les données de niveau d'eau pour la période de trois mois peuvent êtres récupérées pour chaque station.

Les moyennes journalières de niveau d'eau de la période de trois mois seront utilisées pour créer une trame de données pour chaque station à l'aide de Pandas. Pour faciliter la manipulation des données, l'ensemble des trames de données sera rassemblée dans un dictionnaire Python. Ces trames de données seront utilisées pour créer des graphiques interactifs et des tableaux affichant les données historiques de niveau d'eau ainsi qu'une carte interactive présentant les données les plus récentes pour la période temporelle sélectionnée.

Les stations n'ayant aucune donnée de niveau d'eau disponible durant la période choisie seront retirées du jeu de données.

In [6]:
# Récupération des données hydrométriques pour chaque station

# Dictionnaire qui contiendra une trame de données de niveau
# d'eau pour chaque station
hydrometric_data = {}

# Liste des stations n'ayant aucune donnée de niveau d'eau 
stations_without_data = []

# Récupération des données et création des trames de données
for station in stations:

    # Récupération des données de niveaux d'eau
    hydro_data = oafeat.collection_items(
        "hydrometric-daily-mean",
        bbox=bbox,
        datetime=f"{start_date}/{end_date}",
        STATION_NUMBER=station,
    )
    # Création des trames de données s'il y a 
    # des données pour la période choisie
    if hydro_data["features"]:
        # Création d'un dictionnaire dans un format compatible avec Pandas
        historical_data_format = [
            {
                "LATITUDE": el["geometry"]["coordinates"][1],
                "LONGITUDE": el["geometry"]["coordinates"][0],
                **el["properties"],
            }
            for el in hydro_data["features"]
        ]
        # Création d'une trame de données à partir du dictionnaire
        historical_data_df = pd.DataFrame(
            historical_data_format,
            columns=[
                "STATION_NUMBER",
                "STATION_NAME",
                "DATE",
                "LEVEL",
                "LATITUDE",
                "LONGITUDE",
            ],
        )
        historical_data_df = historical_data_df.fillna(value=np.nan)
        # Ajout de la trame de données au dictionnaire hydrometric_data
        if not historical_data_df["LEVEL"].isnull().all():
            # Retrait des lignes à la fin de la trame de données qui ne
            # contiennent pas de données
            while np.isnan(historical_data_df["LEVEL"].iloc[-1]):
                historical_data_df = historical_data_df.drop(
                    historical_data_df.tail(1).index
                )
            # Création d'un index avec la date dans un format datetime
            historical_data_df["DATE"] = pd.to_datetime(
                historical_data_df["DATE"]
            )
            historical_data_df.set_index(["DATE"], inplace=True, drop=True)
            historical_data_df.index = historical_data_df.index.date
            # Ajout de la trame de données au dictionnaire
            hydrometric_data[station] = historical_data_df
        # Si toutes les données sont des NaN, la station sera retirée
        # du jeu de données
        else:
            stations_without_data.append(station)
    # Si aucune donnée n'est disponible pour la période ciblée, la station
    # sera retirée du jeu de données
    else:
        stations_without_data.append(station)

# Retrait des stations sans données de niveau d'eau de la liste des stations
for station in stations_without_data:
    stations.remove(station)

# Générer une erreur si la liste de stations est vide
if not stations:
    raise ValueError(
        "Aucune données de niveau d'eau n'est disponible dans cette période"
        + f" de {num_months} mois pour les stations sélectionnées."
    )

Maintenant que toutes les trames de données sont prêtes, les données peuvent être visualisées sous la forme d'un graphique. Matplotlib peut être utilisé pour générer des graphiques interactifs à l'aide de son API de traitement d'événements.

In [7]:
# Création d'un graphique interactif avec Matplotlib

# Station hydrométrique à afficher sur le graphique
station_displayed_p = stations[1]

# Fonction pour créer un graphique pour la station sélectionnée
def interactive_plot(station):
    # Ajustement de la taille des caractères et de la figure
    params = {
        "legend.fontsize": "14",
        "figure.figsize": (9, 5),
        "axes.labelsize": "14",
        "axes.titlesize": "16",
        "xtick.labelsize": "12",
        "ytick.labelsize": "12",
    }
    plt.rcParams.update(params)
    
    # Création du graphique
    fig, ax = plt.subplots()
    line, = plt.plot(
        hydrometric_data[station].index,
        hydrometric_data[station]["LEVEL"],
        marker="o",
        label="Moyenne journalière",
    )
    plt.legend()
    plt.grid(True, which="both")
    ax.set_title(
        fill(
            "Niveaux d'eau de la station {} ({})".format(
                hydrometric_data[station]["STATION_NAME"][0], station
            ), 60
        )
    )
    ax.set_ylabel("Niveau d'eau (m)")
    ax.set_xlabel("Date")

    # Modification de la graduation et des étiquettes de l'axe des x
    locator = mdates.AutoDateLocator()
    ax.xaxis.set_major_locator(locator)
    plt.gcf().autofmt_xdate()
    
    # Création des annotations à afficher au passage de la souris
    annot = ax.annotate(
        "",
        xy=(0, 0),
        xytext=(-60, -40),
        textcoords="offset points",
        bbox=dict(boxstyle="round", fc="w"),
        arrowprops=dict(arrowstyle="->"),
    )
    annot.set_visible(False)

    return line, annot, ax, fig


# Création de la figure
line, annot, ax, fig = interactive_plot(station_displayed_p)


# Mise à jour des annotations avec les informations du point de donnée
def update_annot(ind):
    # Identification des annotations à afficher
    x, y = line.get_data()
    annot.xy = (x[ind["ind"][0]], y[ind["ind"][0]])
    
    # Ajout de texte à l'annotation (date et niveau d'eau)
    date_x = x[ind["ind"]][0]
    level_y = round(y[ind["ind"]][0], 2)
    text = "{}\nMoyenne journalière: {} m".format(date_x, level_y)
    annot.set_text(text)
    
    # Sélection du niveau de transparence de l'annotation
    annot.get_bbox_patch().set_alpha(0.8)


# Affichage des annotations au passage de la souris
def hover(event):
    vis = annot.get_visible()
    if event.inaxes == ax:
        cont, ind = line.contains(event)
        if cont:
            update_annot(ind)
            annot.set_visible(True)
            fig.canvas.draw_idle()
        else:
            if vis:
                annot.set_visible(False)
                fig.canvas.draw_idle()


# Ajout de la fonction d'affichage des annotations au 
# passage de la souris à la figure               
fig.canvas.mpl_connect("motion_notify_event", hover)

plt.show()