Comment télécharger la vapeur d’eau intégrée depuis ERA5 avec Python

1 kg m-2 ≈ 1 mm d’eau précipitable

tcwv = 25 kg m-2

25 mm = 2.5 cm d’eau précipitable

total_column_water_vapour

tcwv

kg m-2

pip install cdsapi xarray netCDF4 matplotlib cartopy

conda install -c conda-forge cdsapi xarray netcdf4 matplotlib cartopy

~/.cdsapirc

url: https://cds.climate.copernicus.eu/api
key: YOUR_PERSONAL_ACCESS_TOKEN

required licences not accepted

import cdsapi

dataset = "reanalysis-era5-single-levels"

request = {
    "product_type": ["reanalysis"],
    "variable": ["total_column_water_vapour"],
    "year": ["2019"],
    "month": ["08"],
    "day": ["08"],
    "time": [
        "00:00", "01:00", "02:00", "03:00",
        "04:00", "05:00", "06:00", "07:00",
        "08:00", "09:00", "10:00", "11:00",
        "12:00", "13:00", "14:00", "15:00",
        "16:00", "17:00", "18:00", "19:00",
        "20:00", "21:00", "22:00", "23:00",
    ],
    "data_format": "netcdf",
    "download_format": "unarchived",

    # Boîte englobante : Nord, Ouest, Sud, Est
    "area": [50.0, -130.0, 24.0, -65.0],
}

client = cdsapi.Client()
client.retrieve(dataset, request).download("era5_tcwv_conus_20190808.nc")

[North, West, South, East]

[50.0, -130.0, 24.0, -65.0]

import xarray as xr

ds = xr.open_dataset("era5_tcwv_conus_20190808.nc")

print(ds)

tcwv = ds["tcwv"]

tcwv = ds["tcwv"]

print(tcwv)
print(tcwv.attrs)
print(tcwv.coords)

time_name = "valid_time" if "valid_time" in tcwv.coords else "time"

print("Coordonnée temporelle :", time_name)

import matplotlib.pyplot as plt

tcwv0 = tcwv.isel({time_name: 0})

plt.figure(figsize=(10, 6))

tcwv0.plot(
    cmap="viridis",
    cbar_kwargs={"label": "TCWV (kg m$^{-2}$ / mm)"}
)

plt.title(f"Vapeur d’eau totale de colonne ERA5 sur CONUS\n{tcwv0[time_name].values}")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

tcwv0 = tcwv.isel({time_name: 0})

fig = plt.figure(figsize=(11, 6.5))

ax = fig.add_axes(
    [0.07, 0.12, 0.78, 0.76],
    projection=ccrs.PlateCarree()
)

cax = fig.add_axes([0.88, 0.18, 0.025, 0.62])

p = tcwv0.plot(
    ax=ax,
    transform=ccrs.PlateCarree(),
    cmap="viridis",
    add_colorbar=False
)

cbar = fig.colorbar(p, cax=cax)
cbar.set_label("Vapeur d’eau totale de colonne (kg m$^{-2}$ / mm)")

ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linewidth=0.8)
ax.add_feature(cfeature.STATES, linewidth=0.35, edgecolor="white", alpha=0.8)

ax.set_extent([-130, -65, 24, 50], crs=ccrs.PlateCarree())

gl = ax.gridlines(
    draw_labels=True,
    linewidth=0.3,
    alpha=0.4,
    linestyle="--"
)

gl.top_labels = False
gl.right_labels = False

ax.set_title(
    f"Vapeur d’eau totale de colonne ERA5 sur CONUS\n{tcwv0[time_name].values}",
    fontsize=13
)

plt.show()

cax = fig.add_axes([0.88, 0.18, 0.025, 0.62])

[left, bottom, width, height]

cax = fig.add_axes([0.88, 0.22, 0.025, 0.50])

tcwv_daily_mean = tcwv.mean(dim=time_name)

plt.figure(figsize=(10, 6))

tcwv_daily_mean.plot(
    cmap="viridis",
    cbar_kwargs={"label": "TCWV moyenne journalière (kg m$^{-2}$ / mm)"}
)

plt.title("Moyenne journalière ERA5 de la vapeur d’eau totale de colonne sur CONUS\n2019-08-08")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

fig = plt.figure(figsize=(11, 6.5))

ax = fig.add_axes(
    [0.07, 0.12, 0.78, 0.76],
    projection=ccrs.PlateCarree()
)

cax = fig.add_axes([0.88, 0.18, 0.025, 0.62])

p = tcwv_daily_mean.plot(
    ax=ax,
    transform=ccrs.PlateCarree(),
    cmap="viridis",
    add_colorbar=False
)

cbar = fig.colorbar(p, cax=cax)
cbar.set_label("TCWV moyenne journalière (kg m$^{-2}$ / mm)")

ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linewidth=0.8)
ax.add_feature(cfeature.STATES, linewidth=0.35, edgecolor="white", alpha=0.8)

ax.set_extent([-130, -65, 24, 50], crs=ccrs.PlateCarree())

gl = ax.gridlines(
    draw_labels=True,
    linewidth=0.3,
    alpha=0.4,
    linestyle="--"
)

gl.top_labels = False
gl.right_labels = False

ax.set_title(
    "Moyenne journalière ERA5 de la vapeur d’eau totale de colonne sur CONUS\n2019-08-08",
    fontsize=13
)

plt.show()

tcwv_mm = tcwv
tcwv_cm = tcwv / 10.0

tcwv_cm = tcwv / 10.0
tcwv_cm.attrs["long_name"] = "Vapeur d’eau totale de colonne"
tcwv_cm.attrs["units"] = "cm"

tcwv0_cm = tcwv_cm.isel({time_name: 0})

plt.figure(figsize=(10, 6))

tcwv0_cm.plot(
    cmap="viridis",
    cbar_kwargs={"label": "TCWV (cm)"}
)

plt.title(f"Vapeur d’eau totale de colonne ERA5 sur CONUS\n{tcwv0_cm[time_name].values}")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

tcwv_mean = tcwv.mean(dim=["latitude", "longitude"])

plt.figure(figsize=(8, 4))

tcwv_mean.plot(marker="o")

plt.title("Vapeur d’eau totale de colonne ERA5 moyenne sur le domaine")
plt.ylabel("TCWV (kg m$^{-2}$ / mm)")
plt.xlabel("Temps")
plt.grid(True, alpha=0.3)
plt.show()

tcwv_min = tcwv.min(dim=["latitude", "longitude"])
tcwv_mean = tcwv.mean(dim=["latitude", "longitude"])
tcwv_max = tcwv.max(dim=["latitude", "longitude"])

plt.figure(figsize=(8, 4))

plt.plot(tcwv[time_name], tcwv_min, marker="o", label="min")
plt.plot(tcwv[time_name], tcwv_mean, marker="o", label="mean")
plt.plot(tcwv[time_name], tcwv_max, marker="o", label="max")

plt.title("Statistiques de la TCWV ERA5 sur le domaine")
plt.ylabel("TCWV (kg m$^{-2}$ / mm)")
plt.xlabel("Temps")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

lat = 45.8326
lon = -120.5

import pandas as pd

target_time = pd.Timestamp("2019-08-08T20:15:00")

tcwv_at_time = tcwv.sel({time_name: target_time}, method="nearest")

tcwv_point = tcwv_at_time.interp(
    latitude=lat,
    longitude=lon
)

print(float(tcwv_point.values))

import pandas as pd
import xarray as xr

# Exemple de DataFrame
df = pd.DataFrame({
    "latitude": [40.0, 42.5, 45.0],
    "longitude": [-120.0, -115.0, -110.0],
    "acq_date_time": [
        "2019-08-08T18:20:00",
        "2019-08-08T19:45:00",
        "2019-08-08T20:10:00",
    ],
})

df["acq_date_time"] = pd.to_datetime(df["acq_date_time"])

tcwv_values = []

for _, row in df.iterrows():
    tcwv_at_time = tcwv.sel(
        {time_name: row["acq_date_time"]},
        method="nearest"
    )

    value = tcwv_at_time.interp(
        latitude=row["latitude"],
        longitude=row["longitude"]
    )

    tcwv_values.append(float(value.values))

df["tcwv_kg_m2"] = tcwv_values
df["tcwv_cm"] = df["tcwv_kg_m2"] / 10.0

print(df)

import numpy as np

def clear_sky_transmittance_3p9_from_tcwv(
    tcwv_kg_m2,
    view_zenith_deg,
    a=0.04,
    b=0.06
):
    """
    Transmittance claire approximative de premier ordre à 3.9 µm
    à partir de la vapeur d’eau totale de colonne ERA5.

    Parameters
    ----------
    tcwv_kg_m2 : float or array-like
        Vapeur d’eau totale de colonne ERA5 en kg m-2.
        Numériquement équivalente aux mm d’eau précipitable.

    view_zenith_deg : float or array-like
        Angle zénithal de visée du satellite en degrés.

    a, b : float
        Coefficients empiriques. Ils devraient à terme être ajustés
        ou remplacés par une table de correspondance de transfert radiatif.

    Returns
    -------
    T : float or array-like
        Transmittance atmosphérique approximative entre 0 et 1.
    """

    tcwv_cm = tcwv_kg_m2 / 10.0
    mu = np.cos(np.deg2rad(view_zenith_deg))

    tau = a + b * tcwv_cm

    return np.exp(-tau / mu)

view_zenith_deg = 35.0

T_clear = clear_sky_transmittance_3p9_from_tcwv(
    tcwv_kg_m2=tcwv,
    view_zenith_deg=view_zenith_deg
)

T_clear.attrs["long_name"] = "Transmittance claire approximative à 3.9 µm"
T_clear.attrs["units"] = "1"

T0 = T_clear.isel({time_name: 0})

plt.figure(figsize=(10, 6))

T0.plot(
    cmap="viridis",
    vmin=0,
    vmax=1,
    cbar_kwargs={"label": "Transmittance approximative à 3.9 µm"}
)

plt.title(f"Transmittance claire approximative à 3.9 µm\nVZA = {view_zenith_deg}°")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()