Note

Go to the end to download the full example code.

Water FluxΒΆ

This example shows how to calculate horizontal water vapor flux on pressure levels and its vertical integral from the top of the atmosphere to the surface. The input fields are specific humidity, zonal wind, meridional wind, and surface pressure.

Before proceeding with all the steps, first import some necessary libraries and packages

import easyclimate as ecl
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.ticker as mticker

Read sample data

The wind and specific humidity fields are pressure-level data. The surface pressure field is used when integrating the flux through the atmospheric column.

u_data = ecl.open_tutorial_dataset("uwnd_2022_day5")["uwnd"].sortby("lat")
v_data = ecl.open_tutorial_dataset("vwnd_2022_day5")["vwnd"].sortby("lat")
q_data = ecl.open_tutorial_dataset("shum_2022_day5")["shum"].sortby("lat")
slp_data = ecl.open_tutorial_dataset("pres_2022_day5")["pres"].sortby("lat")

uwnd_2022_day5.nc ━━━━━━━━━━━━━━━━━━━━ 100.0% β€’ 2.2/2.2 MB β€’ 24.0 MB/s β€’ 0:00:00
vwnd_2022_day5.nc ━━━━━━━━━━━━━━━━━━━━ 100.0% β€’ 2.2/2.2 MB β€’ 18.0 MB/s β€’ 0:00:00
shum_2022_day5.nc ━━━━━━━━━━━━━━━━━━━━ 100.0% β€’ 1.2/1.2 MB β€’ 14.2 MB/s β€’ 0:00:00
pres_2022_day5.nc ━━━━━━━━━━━━━━━━━ 100.0% β€’ 122.6/122.6 kB β€’ 7.0 MB/s β€’ 0:00:00

Define a small helper for the map background. Keeping this in a function makes the single-level and vertically integrated plots use the same map extent, coastlines, and grid labels.

def draw_base_background(ax, lonlat_range):
    ax.set_extent(lonlat_range, crs=ccrs.PlateCarree())
    ax.coastlines(resolution="50m", linewidth=0.8, color="grey")
    gl = ax.gridlines(
        crs=ccrs.PlateCarree(),
        draw_labels=True,
        linewidth=0.5,
        color="gray",
        alpha=0.4,
        linestyle="--"
    )
    gl.top_labels = False
    gl.right_labels = False
    gl.left_labels = True
    gl.bottom_labels = True
    gl.xlocator = mticker.FixedLocator([80, 100, 120, 140, 160, 180])
    gl.ylocator = mticker.FixedLocator([0, 10, 20, 30, 40, 50, 60, 70])

Horizontal Water Vapor FluxΒΆ

easyclimate.calc_horizontal_water_flux calculates water vapor flux at each pressure level:

\[\frac{1}{g} q \mathbf{V} = \frac{1}{g} (u q\ \mathbf{i} + vq\ \mathbf{j})\]

The result is an xarray.Dataset with two components:

  • qu: zonal water vapor flux

  • qv: meridional water vapor flux

qflux = ecl.calc_horizontal_water_flux(
    specific_humidity_data = q_data,
    u_data = u_data,
    v_data = v_data,
)
qflux

<xarray.Dataset> Size: 3MB
Dimensions:  (time: 5, lon: 144, lat: 73, level: 8)
Coordinates:
  * time     (time) datetime64[ns] 40B 2022-01-01 2022-01-02 ... 2022-01-05
  * lon      (lon) float32 576B 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat      (lat) float32 292B -90.0 -87.5 -85.0 -82.5 ... 82.5 85.0 87.5 90.0
  * level    (level) float32 32B 1e+03 925.0 850.0 700.0 600.0 500.0 400.0 300.0
Data variables:
    qu       (time, level, lat, lon) float32 2MB 0.0004865 ... -1.15e-05
    qv       (time, level, lat, lon) float32 2MB -0.001453 ... -1.976e-05


Plot the horizontal water vapor flux at 850 hPa on 2022-01-03. The vectors show the direction and relative magnitude of water vapor transport.

lonlat_range = [70, 180, 0, 70]

fig, ax = plt.subplots(
    1, 1,
    figsize=(4.5, 5),
    subplot_kw={"projection": ccrs.PlateCarree()},
    constrained_layout=True
)

draw_base_background(ax, lonlat_range)

qflux.isel(time = 2).sel(level = 850).plot.quiver(
    x = "lon", y = "lat", u = "qu", v = "qv",
    transform=ccrs.PlateCarree(),
    regrid_shape = 15,
    scale = 0.1,
    add_guide = False
)

ax.set_title("Water Flux (850hPa, 2022.1.3)")
Water Flux (850hPa, 2022.1.3)
Text(0.5, 1.0, 'Water Flux (850hPa, 2022.1.3)')

Vertically Integrated Water Vapor FluxΒΆ

easyclimate.calc_water_flux_top2surface_integral integrates the horizontal water vapor flux over the atmospheric column:

\[\frac{1}{g} \int_0^{p_s} (q\mathbf{v}),dp\]

The pressure coordinate is given by vertical_dim="level" and vertical_dim_units="hPa". Since the surface pressure field is in Pa, set surface_pressure_data_units="Pa" explicitly.

qflux_integral = ecl.calc_water_flux_top2surface_integral(
    specific_humidity_data = q_data,
    u_data = u_data,
    v_data = v_data,
    surface_pressure_data = slp_data,
    vertical_dim = "level",
    specific_humidity_data_units = "kg/kg",
    surface_pressure_data_units = "Pa",
    vertical_dim_units = "hPa",
    # support backend method selection, e.g., "ncl", "rust", "rust-batch"
    # The following parameters generally do not need to be passed
    method = "rust-block"
)
qflux_integral

<xarray.Dataset> Size: 842kB
Dimensions:  (lon: 144, lat: 73, time: 5)
Coordinates:
  * lon      (lon) float32 576B 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat      (lat) float32 292B -90.0 -87.5 -85.0 -82.5 ... 82.5 85.0 87.5 90.0
  * time     (time) datetime64[ns] 40B 2022-01-01 2022-01-02 ... 2022-01-05
Data variables:
    qu       (time, lat, lon) float64 420kB 6.573 6.839 7.141 ... -1.253 -1.11
    qv       (time, lat, lon) float64 420kB -6.845 -6.551 -6.238 ... 1.855 1.933
Attributes:
    long_name:  Vertical integral of water vapour flux (1/g \int q U dp)
    units:      kg m**-1 s**-1


Plot the vertically integrated flux. Compared with the single-level result, this vector field represents the column-integrated horizontal moisture transport.

lonlat_range = [70, 180, 0, 70]

fig, ax = plt.subplots(
    1, 1,
    figsize=(4.5, 5),
    subplot_kw={"projection": ccrs.PlateCarree()},
    constrained_layout=True
)

draw_base_background(ax, lonlat_range)

qflux_integral.isel(time = 2).plot.quiver(
    x = "lon", y = "lat", u = "qu", v = "qv",
    transform=ccrs.PlateCarree(),
    regrid_shape = 15,
    scale = 5000,
    add_guide = False
)

ax.set_title("Water Flux (Integral, 2022.1.3)")
Water Flux (Integral, 2022.1.3)
Text(0.5, 1.0, 'Water Flux (Integral, 2022.1.3)')

Total running time of the script: (0 minutes 6.828 seconds)