Note

Go to the end to download the full example code.

Geographic Finite Difference#

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

import easyclimate as ecl
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

Then consider obtaining meridional and zonal wind variables in tutorial data

u_data = ecl.tutorial.open_tutorial_dataset("uwnd_202201_mon_mean").sortby("lat").uwnd
v_data = ecl.tutorial.open_tutorial_dataset("vwnd_202201_mon_mean").sortby("lat").vwnd
z_data = ecl.tutorial.open_tutorial_dataset("hgt_202201_mon_mean").sortby("lat").hgt
temp_data = ecl.tutorial.open_tutorial_dataset("air_202201_mon_mean").sortby("lat").air
q_data = ecl.tutorial.open_tutorial_dataset("shum_202201_mon_mean").sortby("lat").shum
msl_data = (
    ecl.tutorial.open_tutorial_dataset("pressfc_202201_mon_mean").sortby("lat").pres
)
pr_data = (
    ecl.tutorial.open_tutorial_dataset("precip_202201_mon_mean").sortby("lat").precip
)

uvdata = xr.Dataset()
uvdata["uwnd"] = u_data
uvdata["vwnd"] = v_data

Obtain data slices on 500hPa isobars for January 2022

uvdata_500_202201 = uvdata.sel(level=500, time="2022-01-01")
z_data_500_202201 = z_data.sel(level=500, time="2022-01-01")
temp_data_500_202201 = temp_data.sel(level=500, time="2022-01-01")

Plotting a sample quiver plot of this data slice

fig, ax = plt.subplots(
    subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)}
)

ax.stock_img()
ax.gridlines(draw_labels=["bottom", "left"], color="grey", alpha=0.5, linestyle="--")
ax.coastlines(edgecolor="black", linewidths=0.5)

uvdata_500_202201.thin(lon=3, lat=3).plot.quiver(
    ax=ax,
    u="uwnd",
    v="vwnd",
    x="lon",
    y="lat",
    # projection on data
    transform=ccrs.PlateCarree(),
)

time = 2022-01-01, level = 500.0 [millibar]

<matplotlib.quiver.Quiver object at 0x7f804dd21150>

Second-order Partial Derivative#

The solution of the second-order partial derivative relies on three functional calculations

easyclimate.calc_lon_laplacian: calculation of the second-order partial derivative term (Laplace term) along longitude.

\[\frac{\partial^2 F}{\partial x^2} = \frac{1}{(R \cos\varphi)^2} \cdot \frac{\partial^2 F}{\partial \lambda^2}\]

uwnd_dlon2 = ecl.calc_lon_laplacian(
    uvdata_500_202201.uwnd, lon_dim="lon", lat_dim="lat"
)

easyclimate.calc_lat_laplacian: calculation of the second-order partial derivative term (Laplace term) along latitude.

\[\frac{\partial^2 F}{\partial y^2} = \frac{1}{R^2} \cdot \frac{\partial^2 F}{\partial \varphi^2}\]

uwnd_dlat2 = ecl.calc_lat_laplacian(uvdata_500_202201.uwnd, lat_dim="lat")

easyclimate.calc_lon_lat_mixed_derivatives: second-order mixed partial derivative terms along longitude and latitude.

\[\frac{\partial^2 F}{\partial x \partial y} = \frac{1}{R^2 \cos\varphi} \cdot \frac{\partial^2 F}{\partial \lambda \partial \varphi}\]

uwnd_dlonlat = ecl.calc_lon_lat_mixed_derivatives(
    uvdata_500_202201.uwnd, lon_dim="lon", lat_dim="lat"
)

Second-order partial derivative term along longitude.

fig, ax = plt.subplots(
    subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)}
)

ax.gridlines(draw_labels=["bottom", "left"], color="grey", alpha=0.5, linestyle="--")
ax.coastlines(edgecolor="black", linewidths=0.5)

uwnd_dlon2.plot.contourf(
    ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={"location": "bottom"}, levels=21
)
ax.set_title("$\\frac{\\partial^2 F}{\\partial x^2}$", fontsize=20)

$\frac{\partial^2 F}{\partial x^2}$

Text(0.5, 1.0326797365031362, '$\\frac{\\partial^2 F}{\\partial x^2}$')

Second-order partial derivative term along latitude.

fig, ax = plt.subplots(
    subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)}
)

ax.gridlines(draw_labels=["bottom", "left"], color="grey", alpha=0.5, linestyle="--")
ax.coastlines(edgecolor="black", linewidths=0.5)

uwnd_dlat2.plot.contourf(
    ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={"location": "bottom"}, levels=21
)
ax.set_title("$\\frac{\\partial^2 F}{\\partial y^2}$", fontsize=20)

$\frac{\partial^2 F}{\partial y^2}$

Text(0.5, 1.048071096550179, '$\\frac{\\partial^2 F}{\\partial y^2}$')

Second-order mixed partial derivative terms along longitude and latitude.

fig, ax = plt.subplots(
    subplot_kw={"projection": ccrs.PlateCarree(central_longitude=180)}
)

ax.gridlines(draw_labels=["bottom", "left"], color="grey", alpha=0.5, linestyle="--")
ax.coastlines(edgecolor="black", linewidths=0.5)

uwnd_dlonlat.plot.contourf(
    ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={"location": "bottom"}, levels=21
)
ax.set_title("$\\frac{\\partial^2 F}{\\partial x \\partial y}$", fontsize=20)

$\frac{\partial^2 F}{\partial x \partial y}$

Text(0.5, 1.02873053875448, '$\\frac{\\partial^2 F}{\\partial x \\partial y}$')

Vorticity and Divergence#

Vorticity and divergence are measures of the degree of atmospheric rotation and volumetric flux per unit volume respectively. For vorticity and divergence in the quasi-geostrophic case, the potential height is used as input data for the calculations. In general, we first calculate the quasi-geostrophic wind.

easyclimate.calc_geostrophic_wind: calculate the geostrophic wind.

\[u_g = - \frac{g}{f} \frac{\partial H}{\partial y}, \ v_g = \frac{g}{f} \frac{\partial H}{\partial x}\]

geostrophic_wind_data_500_202201 = ecl.calc_geostrophic_wind(
    z_data_500_202201, lon_dim="lon", lat_dim="lat"
)

The function easyclimate.calc_vorticity is then used to compute the quasi-geostrophic vorticity.

easyclimate.calc_vorticity: calculate the horizontal relative vorticity term.

\[\zeta = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} + \frac{u}{R} \tan \varphi\]

qg_vor_data_500_202201 = ecl.calc_vorticity(
    u_data=geostrophic_wind_data_500_202201.ug,
    v_data=geostrophic_wind_data_500_202201.vg,
    lon_dim="lon",
    lat_dim="lat",
)

qg_vor_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)

<matplotlib.contour.QuadContourSet object at 0x7f804da15f30>

Similar vorticity for actual winds, but for actual winds rather than quasi-geostrophic winds.

vor_data_500_202201 = ecl.calc_vorticity(
    u_data=uvdata_500_202201["uwnd"],
    v_data=uvdata_500_202201["vwnd"],
    lon_dim="lon",
    lat_dim="lat",
)

vor_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)

<matplotlib.contour.QuadContourSet object at 0x7f804e107010>

In addition, the function easyclimate.calc_divergence calculate the quasi-geostrophic divergence.

\[\mathrm{D} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} - \frac{v}{R} \tan \varphi\]

easyclimate.calc_divergence: calculate the horizontal divergence term.

Quasi-geostrophic divergence

qg_div_data_500_202201 = ecl.calc_divergence(
    u_data=geostrophic_wind_data_500_202201.ug,
    v_data=geostrophic_wind_data_500_202201.vg,
    lon_dim="lon",
    lat_dim="lat",
)

qg_div_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)

<matplotlib.contour.QuadContourSet object at 0x7f804f307370>

Actual divergence

div_data_500_202201 = ecl.calc_divergence(
    u_data=uvdata_500_202201["uwnd"],
    v_data=uvdata_500_202201["vwnd"],
    lon_dim="lon",
    lat_dim="lat",
)

div_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)

<matplotlib.contour.QuadContourSet object at 0x7f804e056350>

Of course, in addition to the built-in finite difference method, the spherical harmonic function mothod can be solved, but you must ensure that it is Global and Regular or Gaussian grid type data.

easyclimate.windspharm.calc_relative_vorticity: calculate the relative vorticity term with the spherical harmonic function mothod.
easyclimate.windspharm.calc_divergence: calculate the horizontal divergence term with the spherical harmonic function mothod.

vor_data_500_202201_windspharm = ecl.windspharm.calc_relative_vorticity(
    u_data=uvdata_500_202201["uwnd"],
    v_data=uvdata_500_202201["vwnd"],
)

vor_data_500_202201_windspharm.sortby("lat").sel(lat=slice(20, 80)).plot.contourf(
    levels=21
)

<matplotlib.contour.QuadContourSet object at 0x7f804dd9b3d0>

div_data_500_202201_windspharm = ecl.windspharm.calc_divergence(
    u_data=uvdata_500_202201["uwnd"],
    v_data=uvdata_500_202201["vwnd"],
)

div_data_500_202201_windspharm.sortby("lat").sel(lat=slice(20, 80)).plot.contourf(
    levels=21
)

<matplotlib.contour.QuadContourSet object at 0x7f804e173250>

Generally speaking, the calculation results of the finite difference method and the spherical harmonic function method are similar. The former does not require global regional data, but the calculation results of the latter are more accurate for high latitude regions.

Advection#

Advection is the process of transport of an atmospheric property solely by the mass motion (velocity field) of the atmosphere; also, the rate of change of the value of the advected property at a given point.

For zonal advection, we can calculate as follows.

\[-u \frac{\partial T}{\partial x}\]

u_advection_500_202201 = ecl.calc_u_advection(
    u_data=uvdata_500_202201["uwnd"], temper_data=temp_data_500_202201
)

u_advection_500_202201.sortby("lat").sel(lat=slice(20, 80)).plot.contourf(levels=21)

<matplotlib.contour.QuadContourSet object at 0x7f804d8df190>

Similarly, the meridional advection can acquire as follows.

\[-v \frac{\partial T}{\partial y}\]

v_advection_500_202201 = ecl.calc_v_advection(
    v_data=uvdata_500_202201["vwnd"], temper_data=temp_data_500_202201
)

v_advection_500_202201.sortby("lat").sel(lat=slice(20, 80)).plot.contourf(levels=21)

<matplotlib.contour.QuadContourSet object at 0x7f804d7d6fe0>

Geographic Finite Difference

Contents

Geographic Finite Difference#

First-order Partial Derivative#

Second-order Partial Derivative#

Vorticity and Divergence#

Advection#

Water Flux#

Water Vapor Flux Divergence#