Geographic Finite Difference

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

Tip

You can download following datasets here:

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

vwnd_202201_mon_mean.nc ━━━━━━━━━━━━━━ 100.0% • 1.2/1.2 MB • 31.8 MB/s • 0:00:00
hgt_202201_mon_mean.nc ━━━━━━━━━━━━ 100.0% • 960.5/960.5   • 29.0 MB/s • 0:00:00
                                             kB
air_202201_mon_mean.nc ━━━━━━━━━━━━━━━ 100.0% • 1.1/1.1 MB • 37.7 MB/s • 0:00:00
shum_202201_mon_mean.nc ━━━━━━━━━━━━ 100.0% • 603.3/603.3  • 32.0 MB/s • 0:00:00
                                              kB
pressfc_202201_mon_mean.nc ━━━━━━━━━━━ 100.0% • 98.7/98.7   • 6.7 MB/s • 0:00:00
                                                kB
precip_202201_mon_mean.nc ━━━━━━━━━━━━ 100.0% • 74.4/74.4   • 6.3 MB/s • 0:00:00
                                                kB

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 0x7fb21a55c510>

First-order Partial Derivative#

Consider the function easyclimate.calc_gradient to compute the gradient of the zonal wind with respect to longitude.

\[\frac{\partial u}{\partial \lambda}\]

The argument dim to the function easyclimate.calc_gradient specifies that the direction of the solution is longitude.

uwnd_dx = ecl.calc_gradient(uvdata_500_202201.uwnd, dim="lon")

uwnd_dx

<xarray.DataArray 'uwnd' (lat: 73, lon: 144)> Size: 42kB
array([[ 0.14999962,  0.13225818,  0.11854815, ...,  0.15927446,
         0.15080655,  0.14274216],
       [ 0.25967765,  0.25645185,  0.24798417, ...,  0.27499998,
         0.27661264,  0.27983665],
       [ 0.33306527,  0.30725718,  0.2846775 , ...,  0.35040367,
         0.34435558,  0.34112835],
       ...,
       [-0.1056456 , -0.12661296, -0.13790333, ..., -0.08024192,
        -0.09354854, -0.11290336],
       [-0.07903278, -0.08870995, -0.09435463, ..., -0.08266127,
        -0.08064544, -0.07741928],
       [-0.0576601 , -0.03749967, -0.02540326, ..., -0.07540298,
        -0.06532216, -0.06693745]], shape=(73, 144), dtype=float32)
Coordinates:
  * lat      (lat) float32 292B -90.0 -87.5 -85.0 -82.5 ... 82.5 85.0 87.5 90.0
  * lon      (lon) float32 576B 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
    time     datetime64[ns] 8B 2022-01-01
    level    float32 4B 500.0


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)

uwnd_dx.plot.contourf(
    ax=ax,
    # projection on data
    transform=ccrs.PlateCarree(),
    # Colorbar is placed at the bottom
    cbar_kwargs={"location": "bottom"},
    levels=21,
)
time = 2022-01-01, level = 500.0 [millibar]
<cartopy.mpl.contour.GeoContourSet object at 0x7fb218e72d90>

Of course, it is also possible to pass in xarray.Dataset directly into the function easyclimate.calc_gradient to iterate through all the variables, so that you can get the gradient of both the zonal and meridional winds with respect to longitude at the same time.

uvwnd_dx = ecl.calc_gradient(uvdata_500_202201, dim="lon")

uvwnd_dx

<xarray.Dataset> Size: 85kB
Dimensions:  (lon: 144, lat: 73)
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     datetime64[ns] 8B 2022-01-01
    level    float32 4B 500.0
Data variables:
    uwnd     (lat, lon) float32 42kB 0.15 0.1323 0.1185 ... -0.06532 -0.06694
    vwnd     (lat, lon) float32 42kB 0.1734 0.1815 0.1895 ... 0.2234 0.2137


However, if one is required to solve for the gradient of the zonal wind with respect to the corresponding distance at each longitude, the function calc_dx_gradient should be used to calculate.

\[\frac{\partial F}{\partial x} = \frac{1}{R \cos\varphi} \cdot \frac{\partial F}{\partial \lambda}\]
uwnd_dlon = ecl.calc_dx_gradient(uvdata_500_202201.uwnd, lon_dim="lon", lat_dim="lat")

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_dlon.plot.contourf(
    ax=ax,
    # projection on data
    transform=ccrs.PlateCarree(),
    # Colorbar is placed at the bottom
    cbar_kwargs={"location": "bottom"},
    levels=21,
)
time = 2022-01-01, level = 500.0 [millibar]
<cartopy.mpl.contour.GeoContourSet object at 0x7fb218f2cd10>

Similarly, use easyclimate.calc_dy_gradient to solve for the gradient of the meridional wind with respect to the corresponding distance at each latitude.

Second-order Partial Derivative#

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

\[\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_dx_laplacian(
    uvdata_500_202201.uwnd, lon_dim="lon", lat_dim="lat"
)
\[\frac{\partial^2 F}{\partial y^2} = \frac{1}{R^2} \cdot \frac{\partial^2 F}{\partial \varphi^2}\]
uwnd_dlat2 = ecl.calc_dy_laplacian(uvdata_500_202201.uwnd, lat_dim="lat")
\[\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_dxdy_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.0480710965501792, '$\\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.

\[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.

\[\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)
time = 2022-01-01, level = 500.0 [millibar]
<matplotlib.contour.QuadContourSet object at 0x7fb219028d90>

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)
time = 2022-01-01, level = 500.0 [millibar]
<matplotlib.contour.QuadContourSet object at 0x7fb2188ba790>

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\]

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)
time = 2022-01-01, level = 500.0 [millibar]
<matplotlib.contour.QuadContourSet object at 0x7fb218781a90>

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)
time = 2022-01-01, level = 500.0 [millibar]
<matplotlib.contour.QuadContourSet object at 0x7fb2186847d0>

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
)
plot geographic finite difference
<matplotlib.contour.QuadContourSet object at 0x7fb2187249d0>

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
)
plot geographic finite difference
<matplotlib.contour.QuadContourSet object at 0x7fb2185d8c50>

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)
time = 2022-01-01, level = 500.0 [millibar]
<matplotlib.contour.QuadContourSet object at 0x7fb2185358d0>

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)
time = 2022-01-01, level = 500.0 [millibar]
<matplotlib.contour.QuadContourSet object at 0x7fb2184da290>

Water Flux#

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

easyclimate.calc_horizontal_water_flux can calculate the horizontal water flux of single layers.

ecl.calc_horizontal_water_flux(
    specific_humidity_data=q_data,
    u_data=uvdata.uwnd,
    v_data=uvdata.vwnd,
)

<xarray.Dataset> Size: 1MB
Dimensions:  (time: 2, lon: 144, lat: 73, level: 8)
Coordinates:
  * time     (time) datetime64[ns] 16B 2022-01-01 2022-02-01
  * 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 673kB 0.3579 0.3982 ... -0.007047
    qv       (time, level, lat, lon) float32 673kB -0.9287 -0.9122 ... 0.001194


The whole layer integral needs to consider the function easyclimate.calc_water_flux_top2surface_integral to calculate.

water_flux_top2surface_integral = ecl.calc_water_flux_top2surface_integral(
    specific_humidity_data=q_data,
    u_data=u_data,
    v_data=v_data,
    surface_pressure_data=msl_data,
    surface_pressure_data_units="millibars",
    specific_humidity_data_units = "kg/kg",
    vertical_dim="level",
    vertical_dim_units="hPa",
)

water_flux_top2surface_integral

<xarray.Dataset> Size: 337kB
Dimensions:  (time: 2, lon: 144, lat: 73)
Coordinates:
  * time     (time) datetime64[ns] 16B 2022-01-01 2022-02-01
  * 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
Data variables:
    qu       (time, lat, lon) float64 168kB 3.565e+03 3.674e+03 ... -2.415e+04
    qv       (time, lat, lon) float64 168kB -2.52e+03 -2.364e+03 ... -1.976e+03
Attributes:
    long_name:  Vertical integral of water vapour flux
    units:      kg m**-1 s**-1


Extracting the entire layer water vapor flux at mid and low latitudes at the 0th time level.

draw_water_flux = (
    water_flux_top2surface_integral.isel(time=0)
    .thin(lon=3, lat=3)
    .sel(lat=slice(-60, 60))
)
draw_pr = pr_data.isel(time=0).sel(lat=slice(-60, 60))

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)

draw_water_flux.plot.quiver(
    ax=ax,
    u="qu",
    v="qv",
    x="lon",
    y="lat",
    transform=ccrs.PlateCarree(),
    zorder=2,
)

draw_pr.plot.contourf(
    ax=ax,
    transform=ccrs.PlateCarree(),
    levels=21,
    cmap="Greens",
    zorder=1,
    cbar_kwargs={"location": "bottom"},
    vmax=20,
)
time = 2022-01-01
<cartopy.mpl.contour.GeoContourSet object at 0x7fb218590d90>

Water Vapor Flux Divergence#

Water vapor flux divergence represents the convergence and divergence of water vapor. There are also two built-in functions to calculate the results of single-layers and whole-layer integration respectively.

\[\nabla \left( \frac{1}{g} q \mathbf{V} \right) = \frac{1}{g} \nabla \cdot \left( q \mathbf{V} \right)\]
divergence_watervaporflux_top2surface_integral = (
    ecl.calc_divergence_watervaporflux_top2surface_integral(
        specific_humidity_data=q_data,
        u_data=u_data,
        v_data=v_data,
        surface_pressure_data=msl_data,
        surface_pressure_data_units="millibars",
        specific_humidity_data_units="grams/kg",
        vertical_dim="level",
        vertical_dim_units="hPa",
    )
)

divergence_watervaporflux_top2surface_integral

<xarray.Dataset> Size: 505kB
Dimensions:  (time: 2, lon: 144, lat: 73)
Coordinates:
  * time     (time) datetime64[ns] 16B 2022-01-01 2022-02-01
  * 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
Data variables:
    qu       (time, lat, lon) float64 168kB 3.565 3.674 3.775 ... -24.05 -24.15
    qv       (time, lat, lon) float64 168kB -2.52 -2.364 ... -3.037 -1.976
    wvdiv    (time, lat, lon) float64 168kB nan nan nan nan ... nan nan nan nan


Extracting the entire layer water vapor flux at mid and low latitudes at the 0th time level.

draw_data = divergence_watervaporflux_top2surface_integral.isel(time=0).sel(
    lat=slice(-60, 60)
).wvdiv

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)

draw_data.plot.contourf(
    ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={"location": "bottom"}, levels=21
)
time = 2022-01-01
<cartopy.mpl.contour.GeoContourSet object at 0x7fb2186a0350>

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

Contents