Last updated on Oct 26, 2024
Demo script for the analyses done in Nakamura and Huang (2018, Science): using xarray
Companion notebook to `demo_script_for_nh2018.ipynb <demo_script_for_nh2018.ipynb>`__, using the xarray interface.
[1]:
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from falwa.xarrayinterface import QGDataset
Load ERA-Interim reanalysis data retrieved from ECMWF server
[2]:
data = xr.open_mfdataset("2005-01-23_to_2005-01-30_[tuv].nc")
Compute LWA diagnostics
The QGDataset automatically recognizes the variable names from ERA data and flips dimensions as required by the underlying QGField objects.
[3]:
qgds = QGDataset(data)
uvtinterp = qgds.interpolate_fields()
refstates = qgds.compute_reference_states()
lwadiags = qgds.compute_lwa_and_barotropic_fluxes()
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
Do scipy interpolation
1826819616 0 4 converged at n = 951
1826819616 0 4 converged at n = 727
1826819616 0 4 converged at n = 950
1826819616 0 4 converged at n = 721
1826819616 0 4 converged at n = 949
1826819616 0 4 converged at n = 716
1826819616 0 4 converged at n = 947
1826819616 0 4 converged at n = 719
1826819616 0 4 converged at n = 947
1826819616 0 4 converged at n = 720
1826819616 0 4 converged at n = 945
1826819616 0 4 converged at n = 721
1826819616 0 4 converged at n = 942
1826819616 0 4 converged at n = 711
1826819616 0 4 converged at n = 944
1826819616 0 4 converged at n = 715
1826819616 0 4 converged at n = 938
1826819616 0 4 converged at n = 703
1826819616 0 4 converged at n = 940
1826819616 0 4 converged at n = 707
1826819616 0 4 converged at n = 940
1826819616 0 4 converged at n = 701
1826819616 0 4 converged at n = 939
1826819616 0 4 converged at n = 707
1826819616 0 4 converged at n = 941
1826819616 0 4 converged at n = 712
1826819616 0 4 converged at n = 940
1826819616 0 4 converged at n = 721
1826819616 0 4 converged at n = 942
1826819616 0 4 converged at n = 718
1826819616 0 4 converged at n = 945
1826819616 0 4 converged at n = 730
1826819616 0 4 converged at n = 943
1826819616 0 4 converged at n = 731
1826819616 0 4 converged at n = 942
1826819616 0 4 converged at n = 734
1826819616 0 4 converged at n = 943
1826819616 0 4 converged at n = 734
1826819616 0 4 converged at n = 941
1826819616 0 4 converged at n = 745
1826819616 0 4 converged at n = 941
1826819616 0 4 converged at n = 744
1826819616 0 4 converged at n = 940
1826819616 0 4 converged at n = 746
1826819616 0 4 converged at n = 941
1826819616 0 4 converged at n = 744
1826819616 0 4 converged at n = 939
1826819616 0 4 converged at n = 749
1826819616 0 4 converged at n = 942
1826819616 0 4 converged at n = 752
1826819616 0 4 converged at n = 939
1826819616 0 4 converged at n = 746
1826819616 0 4 converged at n = 940
1826819616 0 4 converged at n = 734
1826819616 0 4 converged at n = 943
1826819616 0 4 converged at n = 737
1826819616 0 4 converged at n = 944
1826819616 0 4 converged at n = 736
1826819616 0 4 converged at n = 946
1826819616 0 4 converged at n = 738
1826819616 0 4 converged at n = 946
1826819616 0 4 converged at n = 738
1826819616 0 4 converged at n = 944
1826819616 0 4 converged at n = 740
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
line 748: ncforce is None
nd: 61 , jb: 0
nd: 61 , jb: 0
Visualize
[4]:
selected_time = "2005-01-23 00:00:00"
selected_uvtinterp = uvtinterp.sel({ "valid_time": selected_time, "height": 10000 })
selected_refstates = refstates.sel({ "valid_time": selected_time })
selected_lwadiags = lwadiags.sel({ "valid_time": selected_time, "height": 10000 })
3D Variables on one pressure level to display
[5]:
variables_3d = [
(selected_uvtinterp["qgpv"], 'Quasigeostrophic potential vorticity (QGPV)'),
(selected_lwadiags["lwa"], 'Local wave activity (LWA)'),
(selected_uvtinterp["interpolated_u"], 'Interpolated zonal wind (u)'),
(selected_uvtinterp["interpolated_v"], 'Interpolated meridional wind (v)')
]
for variable, name in variables_3d:
plt.figure(figsize=(12, 6))
plt.contourf(variable['xlon'], variable['ylat'], variable, 50, cmap='jet')
if name == 'Local wave activity (LWA)':
plt.axhline(y=0, c='w', lw=30)
plt.colorbar()
plt.ylabel('Latitude (deg)')
plt.xlabel('Longitude (deg)')
plt.title(f"{name} at 10 km | {selected_time}")
Reference states to be displayed on y-z plane
[6]:
variables_yz = [
(selected_refstates["qref"], 'Qref'),
(selected_refstates["uref"], 'Uref'),
(selected_refstates["ptref"], 'PTref')
]
for variable, name in variables_yz:
plt.figure(figsize=(6, 4))
plt.contourf(variable['ylat'], variable['height'], variable, 50, cmap='jet')
plt.axvline(x=0, c='w', lw=2)
plt.xlabel('Latitude (deg)')
plt.ylabel('Pseudoheight (m)')
plt.colorbar()
plt.title(f'{name} | {selected_time}')
Vertically averaged variables to be displayed on x-y plane
[7]:
variables_xy = [
(selected_lwadiags["u_baro"], 'Barotropic zonal wind'),
(selected_lwadiags["lwa_baro"], 'Barotropic LWA'),
(selected_lwadiags["adv_flux_f1"], 'Advective flux F1'),
(selected_lwadiags["adv_flux_f2"], 'Advective flux F2'),
(selected_lwadiags["adv_flux_f3"], 'Advective flux F3'),
(selected_lwadiags["convergence_zonal_advective_flux"], 'Advective flux convergence -Div(F1+F2+F3)'),
(selected_lwadiags["divergence_eddy_momentum_flux"], 'divergence_eddy_momentum_flux'),
(selected_lwadiags["meridional_heat_flux"], 'meridional_heat_flux')
]
for variable, name in variables_xy:
plt.figure(figsize=(12, 6))
plt.contourf(variable['xlon'], variable['ylat'][1:-1], variable[1:-1, :], 50, cmap='jet')
plt.axhline(y=0, c='w', lw=30)
plt.ylabel('Latitude (deg)')
plt.xlabel('Longitude (deg)')
plt.colorbar()
plt.title(f'{name} | {selected_time}')
plt.show()
Write to disk
Write with xarray’s to_netcdf.
[8]:
xr.merge([uvtinterp, refstates, lwadiags]).to_netcdf("2005-01-23_to_2005-01-30_output_xr.nc")
To reduce the file size, compression and/or integer-based value packing can be added by specifying an appropriate encoding for the variables.