A script for creating Figure X of the LOVE19 ESSD Paper¶
- This script assumes it is located inside a directory with the same architecture as the ESSD repository
- The script demonstrates the use of the vertical elements of LOVE19: SODAR-RASS, LIDAR, FlyFOX-V, and FODS on the 12m tower.
In [1]:
# netcdf/numpy/xray/stats
import numpy as np
import xarray as xr
# OS interaction
import os
# import plotting
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import datetime
import pandas as pd
import pywt
import pyfocs
Format plots¶
In [2]:
%matplotlib inline
# Higher resolution figures within the notebook
%config InlineBackend.figure_format = 'retina'
# Set the plot style from the seaborn library
sns.set_style("whitegrid")
sns.set_context("paper", rc={"font.size":8,"axes.titlesize":8,"axes.labelsize":8})
# Define a default color palette (this should be fairly color blind friendly)
flatui = ["#3498db", "#FFBF00", "#95a5a6", "#34495e", "#e74c3c", "#9b59b6",]
sns.set_palette(sns.color_palette(flatui))
Directories¶
The below command should return the path to the LOVE19 ESSD repository example-scripts repository.
In [3]:
dir_notebook = os.getcwd()
dir_essd = os.path.abspath(os.path.relpath(dir_notebook, 'ESSD-repository'))
dir_essd
Out[3]:
In [4]:
# Met data
dir_csat_flux = os.path.join(dir_essd, 'CSAT fluxes')
dir_aws = os.path.join(dir_essd, 'AWS')
# Figures
dir_print = os.path.join(dir_essd, 'figures')
# DTS data (North Simba and Outer Array)
dir_fods_cross = os.path.join(dir_essd, 'FODS', 'FODS-cross')
dir_flyfox_data = os.path.join(dir_essd, 'FODS', 'FlyFOX')
# Ground-based remote sensing
dir_sodar = os.path.join(dir_essd, 'remote-sensing', 'SODAR-RASS')
dir_lidar = os.path.join(dir_essd, 'remote-sensing', 'LIDAR')
Denoising functions¶
FlyFOX-V data are denoised using a hard-threshold wavelet denoising. The hard threshold is determined for a particular critical wavelength.
These rely on the python wavelet package, pywt. This step is not strictly necessary.
In [5]:
# Denoising functions
def lowpassfilter(signal, thresh, wavelet="db8", mode='soft', pad='constant'):
coeff = pywt.wavedec(signal, wavelet, mode=pad)
if mode == 'firm':
coeff[1:] = (pywt.threshold_firm(i, np.min(thresh), np.max(thresh)) for i in coeff[1:])
else:
coeff[1:] = (pywt.threshold(i, value=thresh, mode=mode ) for i in coeff[1:])
reconstructed_signal = pywt.waverec(coeff, wavelet, mode="constant" )
return reconstructed_signal
def lowpassfilter_level(signal, level, wavelet="db8", mode='soft', pad='constant'):
coeff = pywt.wavedec(signal, wavelet, mode=pad)
coeff[-level:] = (np.zeros_like(c) for c in coeff[-level:])
reconstructed_signal = pywt.waverec(coeff, wavelet, mode="constant" )
return reconstructed_signal
def wv_dc(lv, res, omega_0=4):
a_e = 2**lv
length_scale = (a_e * res * np.pi) / (omega_0)
return(length_scale)
def wavelet_filter(ds, waveletname, lvl, dx, level_or_thresh='level', data_var='theta_v'):
# Build the xarray Dataset for the denoised data. We will pass this back to the user.
denoised = xr.Dataset({data_var: (['time', 'z'], np.zeros_like(ds[data_var].values)),},
coords={'time': ds.time,
'z': ds.z.values})
for t in ds.time:
# Construct the signal to be denoised
signal_wnans = ds[data_var].sel(time=t)
signal = signal_wnans.dropna(dim='z')
signal_mean = signal.mean(dim='z')
signal = signal - signal_mean
# Low pass filter using level reduction
if level_or_thresh == 'level':
crit_scale = wv_dc(lvl, dx)
rec = lowpassfilter_level(signal.values, lvl, wavelet=waveletname)
# Low pass filter using thresholding
elif level_or_thresh == 'thresh':
rec = lowpassfilter(signal.values,
lvl,
wavelet=waveletname,
mode='soft')
# Repad with nans that we removed in the first step
if rec.size < signal_wnans.size:
rec = np.append(rec, (np.zeros(signal_wnans.size - rec.size) * np.nan))
# This handles cases where an extra point was added at the end of the data
else:
rec = rec[0:signal_wnans.size]
# Re-add the mean value and assign to the xarray
rec = rec + signal_mean.values
denoised[data_var].loc[dict(time=t)] = rec
if level_or_thresh == 'level':
denoised.attrs['crit_scale'] = crit_scale
return denoised
Dictionary of times¶
In [6]:
# UTC flight times for FlyFOX
flight_times = {'LOVE_Flyfox_190718': ['2019-07-18 3:10', '2019-07-18 6:54'],
'LOVE_Flyfox_190722': ['2019-07-22 3:18', '2019-07-22 6:10'],
'LOVE_Flyfox_190723': ['2019-07-23 2:59', '2019-07-23 6:30'],
'LOVE_Flyfox_190726': ['2019-07-26 3:21', '2019-07-26 6:24'],
}
AWS¶
In [7]:
os.chdir(dir_aws)
aws = xr.open_dataset('AWS_Voi_1min.nc')
DTS Tower data¶
The FODS-cross directory contains the data with the 12m tower. The data for July 18th is gathered as a dask object, pruned to the flight times and loaded into memory, and the tower data selected. Finally, the z-coordinate is updated from a field site relative coordinate to a height above ground.
In [8]:
os.chdir(dir_fods_cross)
ds_ns_hot = xr.open_mfdataset('*2019-07-18*_heated.nc', combine='by_coords')
ds_ns_cold = xr.open_mfdataset('*2019-07-18*_unheated.nc', combine='by_coords')
ds_ns_hot = ds_ns_hot.sel(time=slice(None, '2019-07-18 08:00')).load()
ds_ns_cold = ds_ns_cold.sel(time=slice(None, '2019-07-18 08:00')).load()
In [9]:
ds_tower_cold = ds_ns_cold.where((ds_ns_cold.unheated == 'tower'), drop=True)
ds_tower_hot = ds_ns_hot.where((ds_ns_hot.heated == 'tower'), drop=True)
# Convert the z coordinate to a height above ground.
# Account for the additional 10cm since we didn't measure the exact bottom.
z_tb = 1.448183004 + 0.1
ds_tower_cold['z'] = ds_tower_cold['z'] + z_tb
ds_tower_hot['z'] = ds_tower_hot['z'] + z_tb
In [10]:
# Fiber parameters
rad = 1.34 / 2 * 10**(-3) # radius in meters
volume_per_m = rad ** 2 * np.pi * 1
kg_per_m = 5 / 1000
density = kg_per_m / volume_per_m
density_per_m_fiber = density * volume_per_m
p = 4.5
params = {
'rad': rad,
'crv': 502,
'density': density_per_m_fiber
}
aws_sub = aws.reindex_like(ds_tower_cold.time, method='nearest')
tower_wind_speed = pyfocs.wind_speed.calculate(
ds_tower_hot.cal_temp + 273.15,
ds_tower_cold.cal_temp + 273.15,
p,
(aws_sub['Rlwd (CNR4)'] + aws_sub['Rlwu (CNR4)']),
method='vR20',
params=params)
FlyFOX observations¶
In [11]:
# Open up the July 18th flight
os.chdir(dir_flyfox_data)
ff_18 = xr.open_dataset('flyfox_180719_thetav.nc')
# Denoising parameters
waveletname = 'bior5.5'
level = 5
dx = 0.127
# Denoise using a hard wavelet threshold
denoised_ff18 = wavelet_filter(ff_18, waveletname, level, dx)
crit_scale = denoised_ff18.crit_scale
SODAR¶
In [12]:
os.chdir(dir_sodar)
sodar = xr.open_dataset('LOVE_SODAR_JuneJuly.nc')
LIDAR¶
In [13]:
os.chdir(dir_lidar)
lidar = xr.open_dataset('lidar_stares_83s.nc')
lidar_vad = xr.open_dataset('lidar_vad.nc')
Stacked column figure¶
In [14]:
sns.set_context('paper')
# Prep the figure
fig = plt.figure(figsize=(12, 15))
heights = [1, 1, 1, 1, 3, 1, 1]
widths = [1, 0.025]
nrows = 7
ncols = 2
spec = fig.add_gridspec(ncols=ncols, nrows=nrows,
width_ratios=widths,
height_ratios=heights,
wspace=0.01, hspace=0.16)
ax_maps = []
for row in range(nrows):
ax_maps.append(fig.add_subplot(spec[row, 0]))
ax_cbar_lidar = fig.add_subplot(spec[0, 1])
ax_cbar_flyfox = fig.add_subplot(spec[4, 1])
ax_cbar_sodar_dir = fig.add_subplot(spec[1, 1])
ax_cbar_sodar_sp = fig.add_subplot(spec[2, 1])
ax_cbar_tower_u = fig.add_subplot(spec[6, 1])
# Prep the data
ft = 'LOVE_Flyfox_190718'
tflight1 = pd.Timestamp(flight_times[ft][0])
tflight2 = pd.Timestamp('2019-07-18 04:30')
# tflight2 = pd.Timestamp(flight_times[ft][1])
vmin = denoised_ff18.theta_v.min() - 273.13
vmax = denoised_ff18.theta_v.max() - 273.13
# Prepare the SODAR-RASS data
t1 = tflight1.floor('10min')
t2 = tflight2.ceil('10min')
fl_sodar = sodar.sel(time=slice(t1, t2))
sodar_z_max = 300
sodar_z_min = 0
sodar_z_ticks = np.arange(sodar_z_min, sodar_z_max, 60)
# Prepare the LIDAR data
fl_lidar = lidar.sel(time=slice(t1, t2))
lidar_z_max = 600
lidar_z_min = 0
lidar_z_ticks = np.arange(lidar_z_min, lidar_z_max, 100)
# LIDAR
ax = ax_maps[0]
ax.set_facecolor('0.9')
color = lidar.mean_w.T.values
im = ax.pcolormesh(lidar.time, lidar.z, color,
cmap=plt.get_cmap('RdBu_r'),
linewidth=0,
rasterized=True,
vmin=-1, vmax=1)
im.set_edgecolor('face')
ax.set_ylabel('Height (m)')
ax.set_ylim(lidar_z_min, lidar_z_max)
ax.set_yticks(lidar_z_ticks)
ax.set_xlim(pd.Timestamp(tflight1), pd.Timestamp(tflight2))
cbar = fig.colorbar(im, cax=ax_cbar_lidar)
cbar.ax.set_ylabel('$w~ (ms^{-1}$)', rotation=270, labelpad=25)
ax.set_title('a) LIDAR Vertical Wind', loc='left')
# Direction
ax = ax_maps[1]
ax.set_facecolor('0.9')
color = fl_sodar.wind_direction.sel(z=slice(sodar_z_min, sodar_z_max))
im = ax.pcolormesh(color.time, color.z, color.T,
cmap=plt.get_cmap('twilight'),
linewidth=0,
rasterized=True,
vmin=0, vmax=360)
im.set_edgecolor('face')
ax.set_ylabel('Height (m)')
ax.set_ylim(sodar_z_min, sodar_z_max)
ax.set_yticks(sodar_z_ticks)
ax.set_xlim(tflight1, tflight2)
cbar = fig.colorbar(im, cax=ax_cbar_sodar_dir)
cbar.ax.set_ylabel('Direction [-]', rotation=270, labelpad=25)
cbar.set_ticks([0, 90, 180, 270, 360])
cbar.ax.set_yticklabels(['N', 'E', 'S', 'W', 'N'])
ax.set_title('b) SODAR Wind Direction', loc='left')
# Speed
ax = ax_maps[2]
ax.set_facecolor('0.9')
color = fl_sodar.wind_speed.sel(z=slice(sodar_z_min, sodar_z_max))
color = color.where(color < 10)
im = ax.pcolormesh(color.time, color.z, color.T.values,
cmap=plt.get_cmap('cividis'), vmin=0, vmax=5,
linewidth=0,
rasterized=True)
ax.plot([t1, t2], [40, 40], 'k', linewidth=2)
ax.set_ylabel('Height (m)')
ax.set_ylim(sodar_z_min, sodar_z_max)
ax.set_yticks(sodar_z_ticks)
cbar = fig.colorbar(im, cax=ax_cbar_sodar_sp)
cbar.ax.set_ylabel('U [ms$^{-1}$]', rotation=270, labelpad=25)
ax.set_title('c) SODAR Wind Speed', loc='left')
# Plot the SODAR temperatures
ax = ax_maps[3]
ax.set_facecolor('0.9')
color = fl_sodar.temperature.sel(z=slice(sodar_z_min, sodar_z_max))
color = color.where(color > 0)
ax.pcolormesh(color.time, color.z, color.T.values,
cmap=plt.get_cmap('viridis'), vmin=vmin, vmax=vmax,
linewidth=0,
rasterized=True)
ax.set_ylabel('Height (m)')
ax.set_ylim(sodar_z_min, sodar_z_max)
ax.set_yticks(sodar_z_ticks)
ax.set_xlim(tflight1, tflight2)
ax.set_title('d) SODAR-RASS Temperature', loc='left')
# Plot the denoised FlyFOX data
ax = ax_maps[4]
im = ax.pcolormesh(denoised_ff18.time, denoised_ff18.z,
denoised_ff18.theta_v.rolling(time=3, center=True).mean().T - 273.13,
cmap=plt.get_cmap('viridis'), vmin=vmin, vmax=vmax,
linewidth=0,
rasterized=True,)
ax.set_ylabel('Height (m)')
ax.set_xlabel('Time (UTC = local time - 2)')
fig.autofmt_xdate()
ax.set_ylim(0, 210)
ax.set_yticks(np.arange(0, 220, 20))
ax.set_title('e) FlyFOX', loc='left')
ax.set_xlim(tflight1, tflight2)
plt.colorbar(im, cax=ax_cbar_flyfox, extend='both')
ax_cbar_flyfox.set_ylabel(r'$\theta$ [$^{\circ}C$]', rotation=270, labelpad=25)
# Plot the smoothed DTS tower data
ax = ax_maps[5]
tower_smooth = ds_tower_cold.sel(time=slice(tflight1, tflight2))
im = ax.pcolormesh(tower_smooth.time,
tower_smooth.z + z_tb,
(tower_smooth.cal_temp + (tower_smooth.z + z_tb) * 0.01).T.values,
cmap=plt.get_cmap('viridis'), vmin=vmin, vmax=vmax,
linewidth=0,
rasterized=True,)
ax.set_xlim(tflight1, tflight2)
ax.set_ylim(ds_tower_cold.z.min() + z_tb, ds_tower_cold.z.max() + z_tb)
ax.set_ylabel('Height [m]')
ax.set_title('f) Tower FODS - Temperature', loc='left')
# Plot the FODS of wind speed
ax = ax_maps[6]
tower_smooth = tower_wind_speed.sel(time=slice(tflight1, tflight2))
im = ax.pcolormesh(tower_smooth.time,
tower_smooth.z + z_tb,
tower_smooth.T.values,
cmap=plt.get_cmap('Oranges'), vmin=0., vmax=2.5,
linewidth=0,
rasterized=True,)
ax.set_xlim(tflight1, tflight2)
ax.set_ylim(ds_tower_cold.z.min() + z_tb, ds_tower_cold.z.max() + z_tb)
ax.set_ylabel('Height [m]')
ax.set_title('g) Tower FODS - Wind Speed', loc='left')
plt.colorbar(im, cax=ax_cbar_tower_u, extend='max')
ax_cbar_tower_u.set_ylabel('FODS Wind Speed [$ms^{-1}$]', rotation=270, labelpad=25)
os.chdir(dir_print)
fig.savefig('LOVE-Flyfox-190718_column-example.morning-only.pdf', dpi=300, bbox_inches='tight')
fig.savefig('LOVE-Flyfox-190718-column-example.morning-only.png', dpi=300, bbox_inches='tight')
Some quick example plots of the VAD data from the LIDAR¶
In [15]:
sns.set_context('paper')
# Prep the figure
fig = plt.figure(figsize=(12, 4))
widths = [1, 0.025]
nrows = 1
ncols = 2
spec = fig.add_gridspec(ncols=ncols, nrows=nrows,
width_ratios=widths,
wspace=0.01, hspace=0.16)
ax_vad = fig.add_subplot(spec[0, 0])
ax_cbar = fig.add_subplot(spec[0, 1])
# Prep the data
ft = 'LOVE_Flyfox_190718'
tflight1 = pd.Timestamp(flight_times[ft][0])
tflight2 = pd.Timestamp('2019-07-18 04:30')
# Prepare the LIDAR data
fl_lidar = lidar_vad.sel(time=slice(t1, t2))
lidar_z_max = 600
lidar_z_min = 0
lidar_z_ticks = np.arange(lidar_z_min, lidar_z_max, 100)
# LIDAR
ax = ax_vad
ax.set_facecolor('0.9')
color = lidar_vad.phi.T.values
im = ax.pcolormesh(lidar_vad.time, lidar_vad.z, color,
cmap=plt.get_cmap('twilight'),
linewidth=0,
rasterized=True,
vmin=0, vmax=360)
im.set_edgecolor('face')
ax.set_ylabel('Height (m)')
ax.set_ylim(lidar_z_min, lidar_z_max)
ax.set_yticks(lidar_z_ticks)
ax.set_xlim(pd.Timestamp(tflight1), pd.Timestamp(tflight2))
cbar = fig.colorbar(im, cax=ax_cbar)
cbar.ax.set_ylabel('Direction (-)', rotation=270, labelpad=25)
cbar.set_ticks([0, 90, 180, 270, 360])
ax.set_title('LIDAR VAD $\phi$', loc='left')
Out[15]:
In [16]:
sns.set_context('paper')
# Prep the figure
fig = plt.figure(figsize=(12, 4))
widths = [1, 0.025]
nrows = 1
ncols = 2
spec = fig.add_gridspec(ncols=ncols, nrows=nrows,
width_ratios=widths,
wspace=0.01, hspace=0.16)
ax_vad = fig.add_subplot(spec[0, 0])
ax_cbar = fig.add_subplot(spec[0, 1])
# Prep the data
ft = 'LOVE_Flyfox_190718'
tflight1 = pd.Timestamp(flight_times[ft][0])
tflight2 = pd.Timestamp('2019-07-18 04:30')
# Prepare the LIDAR data
fl_lidar = lidar_vad.sel(time=slice(t1, t2))
lidar_z_max = 600
lidar_z_min = 0
lidar_z_ticks = np.arange(lidar_z_min, lidar_z_max, 100)
# LIDAR
ax = ax_vad
ax.set_facecolor('0.9')
color = lidar_vad.mean_U.T.values
im = ax.pcolormesh(lidar_vad.time, lidar_vad.z, color,
cmap=plt.get_cmap('cividis'),
linewidth=0,
rasterized=True,
vmin=0, vmax=5)
im.set_edgecolor('face')
ax.set_ylabel('Height (m)')
ax.set_ylim(lidar_z_min, lidar_z_max)
ax.set_yticks(lidar_z_ticks)
ax.set_xlim(pd.Timestamp(tflight1), pd.Timestamp(tflight2))
cbar = fig.colorbar(im, cax=ax_cbar)
cbar.ax.set_ylabel('Direction (-)', rotation=270, labelpad=25)
cbar.set_ticks([0, 90, 180, 270, 360])
ax.set_title('LIDAR VAD $\phi$', loc='left')
Out[16]:
CSAT data¶
In [18]:
os.chdir(dir_csat_flux)
csat_1min = xr.open_dataset('LOVE19_CSAT_1min_fluxes.v2020-12-03.nc')
In [19]:
# Prep the data
ft = 'LOVE_Flyfox_190718'
tflight1 = pd.Timestamp(flight_times[ft][0])
tflight2 = pd.Timestamp('2019-07-18 04:30')
ds_csat = csat_1min.sel(time=slice(tflight1, tflight2))
fig, ax = plt.subplots(1, 1)
for n in ds_csat.names:
ax.plot(ds_csat.time, ds_csat['u*'].sel(names=n.names), '-o', label=n.values)
ax.legend()
ax.set_ylabel('$u* (m~ s^{-1})$')
ax.set_ylim(0, 0.25)
fig.autofmt_xdate()
fig, ax = plt.subplots(1, 1)
for n in ds_csat.names:
ax.plot(ds_csat.time, ds_csat['U_scalar_mean_unrot'].sel(names=n.names), '-o', label=n.values)
ax.legend()
ax.set_ylabel('$U (m~ s^{-1})$')
ax.set_ylim(0, 3)
fig.autofmt_xdate()
fig, ax = plt.subplots(1, 1)
for n in ds_csat.names:
ax.plot(ds_csat.time, ds_csat['phi'].sel(names=n.names), '-o', label=n.values)
ax.legend()
ax.set_ylabel('$\phi~ (-)$')
fig.autofmt_xdate()
ax.set_ylim(0, 360)
fig, ax = plt.subplots(1, 1)
for n in ds_csat.names:
ax.plot(ds_csat.time, ds_csat['H'].sel(names=n.names), '-o', label=n.values)
ax.set_ylabel('$Q_h~ (Wm^{-2})$')
ax.legend()
ax.set_ylim(-30, 10)
fig.autofmt_xdate()
Generally there is a jumpiness that is probably not real. I suspect that the fog/condensation plays a role.
In [ ]: