# Libraries required for this tutorial...

from datetime import datetime, timedelta
import numpy as np
import pandas as pd
import os
import xarray as xr
import xwrf
import glob

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.gridspec as gridspec

from matplotlib.ticker import MultipleLocator, FormatStrFormatter  # ticker spacing
from matplotlib.colors import Normalize
from mpl_toolkits.axes_grid1 import make_axes_locatable

import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

import metpy
import metpy.calc as mpcalc
from metpy.plots import colortables, Hodograph, SkewT
from metpy.units import units
import hvplot.xarray
import holoviews as hv

import numpy as np
import metpy.calc as mcalc
from metpy.units import units

import ipywidgets as ipw
import panel as pn
import pandas as pd
import panel.widgets as pnw

import warnings
warnings.filterwarnings('ignore')

hv.extension("bokeh")

import panel as pn
pn.extension()

ERROR 1: PROJ: proj_create_from_database: Open of /opt/conda/share/proj failed

Divergence/Convergence

Next, let’s explore the influence of topography by examining divergence within the LASSO simulation

# Set path to LASSO data
path_staging = "/data/project/ARM_Summer_School_2024_Data/lasso_tutorial/cacti/lasso-cacti"  # path on Jupyter
file_list = sorted(glob.glob(f'{path_staging}/20190129/eda09/base/les/subset_d3/corlasso_met_*'))

# Start DASK
from distributed import Client

client = Client('tcp://127.0.0.1:46015')
client

Client

Client-a1b1ef61-19dd-11ef-8333-ee815c04dfcb

Connection method: Direct
Dashboard: http://127.0.0.1:8787/status

Scheduler Info

Scheduler

Scheduler-3ab87bbf-1fe6-45b7-9e92-e8f22460d08f

Comm: tcp://127.0.0.1:46015 Workers: 4
Dashboard: http://127.0.0.1:8787/status Total threads: 8
Started: 3 minutes ago Total memory: 755.55 GiB

Workers

Worker: 0
Comm: tcp://127.0.0.1:37853 Total threads: 2
Dashboard: http://127.0.0.1:45763/status Memory: 188.89 GiB
Nanny: tcp://127.0.0.1:39841
Local directory: /tmp/dask-scratch-space/worker-3n1js6k3
Tasks executing: Tasks in memory:
Tasks ready: Tasks in flight:
CPU usage: 2.0% Last seen: Just now
Memory usage: 652.61 MiB Spilled bytes: 0 B
Read bytes: 7.59 kiB Write bytes: 7.59 kiB

Worker: 1
Comm: tcp://127.0.0.1:35381 Total threads: 2
Dashboard: http://127.0.0.1:42311/status Memory: 188.89 GiB
Nanny: tcp://127.0.0.1:36569
Local directory: /tmp/dask-scratch-space/worker-qjlze_x0
Tasks executing: Tasks in memory:
Tasks ready: Tasks in flight:
CPU usage: 2.0% Last seen: Just now
Memory usage: 507.37 MiB Spilled bytes: 0 B
Read bytes: 7.24 kiB Write bytes: 7.24 kiB

Worker: 2
Comm: tcp://127.0.0.1:33441 Total threads: 2
Dashboard: http://127.0.0.1:38823/status Memory: 188.89 GiB
Nanny: tcp://127.0.0.1:40699
Local directory: /tmp/dask-scratch-space/worker-4x7tpqr8
Tasks executing: Tasks in memory:
Tasks ready: Tasks in flight:
CPU usage: 2.0% Last seen: Just now
Memory usage: 657.16 MiB Spilled bytes: 0 B
Read bytes: 8.88 kiB Write bytes: 8.88 kiB

Worker: 3
Comm: tcp://127.0.0.1:40993 Total threads: 2
Dashboard: http://127.0.0.1:40543/status Memory: 188.89 GiB
Nanny: tcp://127.0.0.1:36477
Local directory: /tmp/dask-scratch-space/worker-irzhla9u
Tasks executing: Tasks in memory:
Tasks ready: Tasks in flight:
CPU usage: 0.0% Last seen: Just now
Memory usage: 497.05 MiB Spilled bytes: 0 B
Read bytes: 10.17 kiB Write bytes: 10.17 kiB 
# Open all 73 LASSO time steps
ds = xr.open_mfdataset(file_list)
ds

<xarray.Dataset> Size: 430GB
Dimensions:               (Time: 73, south_north: 865, west_east: 750,
                           bottom_top: 149)
Coordinates:
  * Time                  (Time) datetime64[ns] 584B 2019-01-29T06:00:00 ... ...
    XLONG                 (south_north, west_east) float32 3MB dask.array<chunksize=(865, 750), meta=np.ndarray>
    XLAT                  (south_north, west_east) float32 3MB dask.array<chunksize=(865, 750), meta=np.ndarray>
    XTIME                 (Time) float32 292B dask.array<chunksize=(1,), meta=np.ndarray>
Dimensions without coordinates: south_north, west_east, bottom_top
Data variables: (12/52)
    ITIMESTEP             (Time) int32 292B dask.array<chunksize=(1,), meta=np.ndarray>
    MUTOT                 (Time, south_north, west_east) float32 189MB dask.array<chunksize=(1, 865, 750), meta=np.ndarray>
    HGT                   (Time, south_north, west_east) float32 189MB dask.array<chunksize=(1, 865, 750), meta=np.ndarray>
    HAMSL                 (Time, bottom_top, south_north, west_east) float32 28GB dask.array<chunksize=(1, 50, 289, 250), meta=np.ndarray>
    P_HYD                 (Time, bottom_top, south_north, west_east) float32 28GB dask.array<chunksize=(1, 50, 289, 250), meta=np.ndarray>
    PRESSURE              (Time, bottom_top, south_north, west_east) float32 28GB dask.array<chunksize=(1, 50, 289, 250), meta=np.ndarray>
    ...                    ...
    MULFC                 (Time, south_north, west_east) float32 189MB dask.array<chunksize=(1, 865, 750), meta=np.ndarray>
    MULNB                 (Time, south_north, west_east) float32 189MB dask.array<chunksize=(1, 865, 750), meta=np.ndarray>
    MULPL                 (Time, south_north, west_east) float32 189MB dask.array<chunksize=(1, 865, 750), meta=np.ndarray>
    MUCAPE                (Time, south_north, west_east) float32 189MB dask.array<chunksize=(1, 865, 750), meta=np.ndarray>
    MUCIN                 (Time, south_north, west_east) float32 189MB dask.array<chunksize=(1, 865, 750), meta=np.ndarray>
    REFL_10CM_MAX         (Time, south_north, west_east) float32 189MB dask.array<chunksize=(1, 865, 750), meta=np.ndarray>
Attributes: (12/39)
    DX:                          500.0
    DY:                          500.0
    SIMULATION_START_DATE:       2019-01-29_06:00:00
    WEST-EAST_GRID_DIMENSION:    751
    SOUTH-NORTH_GRID_DIMENSION:  866
    BOTTOM-TOP_GRID_DIMENSION:   150
    ...                          ...
    doi_isPartOf_lasso-cacti:    https://doi.org/10.5439/1905789
    doi_isDocumentedBy:          https://doi.org/10.2172/1905845
    doi_thisFileType:            https://doi.org/10.5439/1905819
    history:                     processed by user d3m088 on machine cirrus89...
    filename_user:               corlasso_met_2019012900eda09d3_base_M1.m1.20...
    filename_storage:            corlassomet2019012900eda09d3baseM1.m1.201901...
# Grab 10 meter wind data
u_all = ds.UMET[:,10,:,:]
v_all = ds.VMET[:,10,:,:]

# Set the grid spacing of the model data (100 m) 
dx = 100 * units('m')
dy = 100 * units('m')

# Create a divergence array
divergence = np.zeros((u_all.shape[0], u_all.shape[1], u_all.shape[2]))

# Loop through all files and calculate divergence at each time step
for t in np.arange(0, u_all.shape[0], 1): 
    # Wind at each time step
    u = u_all[t,:,:] * units('m/s')
    v = v_all[t,:,:] * units('m/s')

    # Calculate divergence
    div = mcalc.divergence(u,v,dx=dx,dy=dy)
    divergence[t,:,:] = div.values

# Add divergence to the data array
ds["Divergence"]=(['Time', 'south_north', 'west_east'],  divergence)

# Subset the divergence data and compute
ds_subset = ds[["Divergence"]].compute()

# Create plot object for terrain data
terrain_plot = ds.HGT.isel(Time=0).hvplot.contour(x='XLONG',
                                     y='XLAT',
                                     levels=10,
                                     geo=True,
                                     rasterize=True,
                                     title="Terrain Height",
                                     widget_location="bottom",
                                                 height=400,
                                                 width=500,)

# Create plot object for divergence data

divergence_plot = ds_subset.Divergence.hvplot.quadmesh(x='XLONG', y='XLAT', groupby='Time', clim=(-0.01, 0.01),  # sets colormap limits
    widget_type="scrubber", cmap='Coolwarm',
                              geo=True,
    widget_location="bottom", rasterize=True, height=400, width=500)

# Create interactive plot with time slider

#(terrain_plot + divergence_plot).cols(1)
divergence_plot

We can see two features of interest: a band of convergence that is consistently present and strong convergence along what appear to be cold pools associated with convection. Now, let’s plot a time step to explore each.

# Plot a single time step

cf = plt.contourf(ds.XLONG, ds.XLAT, ds['Divergence'].isel(Time=44), cmap='seismic', levels=np.linspace(-0.05,0.05,15))
plt.contour(ds.XLONG, ds.XLAT, ds.HGT[0,:,:], levels=[500, 1000, 1500], colors='k')
#plt.barbs(ds.XLONG[::40, ::40], ds.XLAT[::40, ::40], u[::40, ::40], v[::40, ::40])
plt.xlim(-65.5, -64.5)
plt.ylim(-33, -32)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Divergence and Terrain')
plt.colorbar(cf, label='Divergence (1/s)')

<matplotlib.colorbar.Colorbar at 0x7f1e41206290>
../_images/f3cb1218def8e3f9272f352fa2f885fa933940252f9a0d72523b4b9b30ac3c39.png 
# Plot a single time step

cf = plt.contourf(ds.XLONG, ds.XLAT, ds['Divergence'].isel(Time=54), cmap='seismic', levels=np.linspace(-0.05,0.05,15))
plt.contour(ds.XLONG, ds.XLAT, ds.REFL_10CM_MAX[54,:,:], levels=[50, 60], colors='green')
#plt.barbs(ds.XLONG[::40, ::40], ds.XLAT[::40, ::40], u[::40, ::40], v[::40, ::40])
plt.xlim(-65.2, -64.2)
plt.ylim(-33, -32)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Divergence and Reflectivity')
plt.colorbar(cf, label='Divergence (1/s)')

<matplotlib.colorbar.Colorbar at 0x7f1e4105f090>
../_images/e56ee945870bd88e41792db0b31826d245fa79f0b6ae159d963cc06f67622c8e.png

So it appears that the terrain is playing a role in initiating convection with convergence along the ridge of the mountain range adding lift to initiate convection. Subsequently, outflow from the initial storms generates more convergence and lift to help initiate more convection.

ds.WA[44,10,:,:].max().compute()

<xarray.DataArray 'WA' ()> Size: 4B
array(4.2769575, dtype=float32)
Coordinates:
    Time     datetime64[ns] 8B 2019-01-29T17:00:00
    XTIME    float32 4B 660.0