Supplementary material 6/6

This notebook serves as supplementary material to the manuscript on low-level jets submitted to Wind Energy Science. The supplementary material consists of 6 notebooks:

  1. A very quick overview of the measurement data from the LiDARs.
  2. Overview of the ERA5 data and especially the procedure of aligning it with the observations.
  3. Some error statistics on wind speed from ERA5 and observations
  4. Main results on low-level jets
  5. Reconstruction of seasonal cycle using a machine learning package
  6. Spatial figures from ERA5 (data-intensive).

In this notebook (6/6), I create some of the spatial plots that are based on large datasets of ERA5 data.

In [1]:
import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
from cartopy import crs as ccrs
from cartopy.feature import NaturalEarthFeature as NEF

# Low-level jet detection
import sys
sys.path.append('../../code/utils')
import lljs

plt.xkcd()
Out[1]:
<matplotlib.rc_context at 0x2aec26457e80>
In [2]:
# Read in coarse grid SLP data
msl = xr.open_mfdataset('../../data/external/syn_class/era*')

# Read in high resolution wind speed data for LLJ detection; Reverse order of levels!
era5 = xr.open_mfdataset('../../data/external/ECMWF/ERA5/era5_*_ml.nc')[['u', 'v']].sel(level=slice(None,124,-1))
wspd = (era5.u**2+era5.v**2)**.5
falloff = lljs.detect_llj_vectorized(wspd.values, axis=1, output='falloff', mask_inv=True, inverse=False)
lljs = xr.DataArray((falloff>2).astype(int), 
                    coords={'time': wspd.time, 'latitude': wspd.latitude, 'longitude': wspd.longitude}, 
                    dims=['time', 'latitude', 'longitude'])
    
# Read in lamb weather type data (I made this in another project)
lwt = pd.read_csv('../../data/external/syn_class/labels_LWT_Jones.csv', index_col=0, parse_dates=True).LWT

# LWTs are automatically sorted on groupby, I want to sort them in a different order:
mapping = ['CN', 'N', 'AN',
           'CNE', 'NE', 'ANE',
           'CE', 'E', 'AE',
           'CSE', 'SE', 'ASE',
           'CS', 'S', 'AS',
           'CSW', 'SW', 'ASW',
           'CW', 'W', 'AW',
           'CNW', 'NW', 'ANW',
           'C', 'U', 'A']
In [3]:
def geowind(mslp, lat=55, dx=2.5, dy=2.5):
    """ Estimate geostrophic wind field based on mean sea level pressure. """
    # Note that I'm using mslp rather than msl to distinguish a single time step (2D) from the full dataset (3D)

    f = 1.e-4 # coriolis parameter (approximate value)
    rho = 1.23 # density (approximate value)
    r = 6371000 # radius of earth
    dx = dx*np.pi/180.*r*np.cos(lat*np.pi/180.) # convert 2.5 degrees to meters at 55N
    dy = -dy*np.pi/180.*r # latitude decreases along axis
    Ug = -1/(f*rho)*np.gradient(mslp.values, dy, axis=0) # np.gradient only works on numpy arrays
    Vg = 1/(f*rho)*np.gradient(mslp.values, dx, axis=1) # .values extracts the bare numpy array from xarray
    return Ug, Vg
In [4]:
def add_windturbine(fig,ax):
    """ This function adds an overlay axes and plots a wind turbine on it. """
    img=plt.imread('../../reports/figures/wind_turbine_black.png')
    aspect = img.shape

    # Create new axes
    pos1 = ax.get_position()
    h = pos1.height/2
    w = h*aspect[1]/aspect[0]
    a = pos1.x1-0.5*w # x-origin
    b = pos1.y0 # y-origin
    ax2 = fig.add_axes([a,b,w,h])

    # Display the turbine
    ax2.imshow(img)
    ax2.set_axis_off()
    return ax2
In [15]:
# First plot one overall average figure

# monkey patch to fix incompatibility between cartopy and matplotlib 
# https://github.com/SciTools/cartopy/issues/1120#issuecomment-424418760
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh

# Open figure
fig = plt.figure(figsize=(10, 4), dpi=200)
ax0 = fig.add_axes([0.05, 0.05, 0.3, 0.9])
ax = fig.add_axes([0.4, 0.05, 0.5, 0.9], projection=ccrs.LambertConformal(3.43, 52.85, standard_parallels=(45,55)))

#################################################################
# Plot illustration on first subplot
def vwind(z, ug=10, vg=0, f=1e-4, K=.2):
    """ v-wind component of Ekman profiles"""   
    gamma = np.sqrt(f/(2*K))
    return vg - vg*np.exp(-gamma*z)*np.cos(gamma*z) + ug*np.exp(-gamma*z)*np.sin(gamma*z)

def logprofile(z, ust=0.3, z0=0.01):
    """ Logarithmic wind profile """
    return ust/0.4*np.log(z/z0)

# Compute normal (logprofile) and low-level jet (Ekman v-profile) wind profiles
z = np.arange(200)+1.
v_normal = logprofile(z=z,ust=0.25, z0=0.7)
v_lljet = vwind(z=z, ug=7, vg=2, K=.1)


# ax0 = fig.add_subplot(121)
ax0.plot(v_normal, z, label="'standard' profile")
ax0.plot(v_lljet, z, label='Low-level jet')
ax0.set_xlim(-1, 5)
ax0.set_ylim(0, 220)
ax0.plot([2, 3.7], [50, 50],'k')
ax0.plot([2, 2], [46, 54],'k')
ax0.plot([3.75, 3.75], [46, 54],'k')
ax0.text(2.35, 57, 'falloff')
ax0.legend(loc=2, fancybox=True, framealpha=0.3)
ax0.set_xticks([])
ax0.set_yticks([])
ax0.set_xlabel('wind speed')
ax0.set_ylabel('height')
ax0.spines['right'].set_visible(False)
ax0.spines['top'].set_visible(False)
ax2 = add_windturbine(fig, ax0)


#################################################################
# Plot map of overall LLJ occurrence on second subplot
# ax = fig.add_subplot(122, projection=ccrs.LambertConformal(3.43, 52.85, standard_parallels=(45,55)))
ax.set_extent([-7, 15, 50, 60])
ax.outline_patch.set_visible(False)

# Add country borders
borders = NEF(category='cultural', name='admin_0_countries', scale='10m', facecolor='none') # 110, 50 or 10
ax.add_feature(borders, edgecolor='gray')

# Plot LLJ spatial distribution
mean_lljs = lljs.mean(dim='time')*100
lon = lljs.longitude
lat = lljs.latitude
mesh = ax.pcolormesh(lon, lat, mean_lljs.values, transform=ccrs.PlateCarree(), cmap='GnBu')

# Colobar
cbar_ax = fig.add_axes([0.84, 0.1, 0.03, 0.8])
cbar = fig.colorbar(mesh, cax=cbar_ax, label='low-level jet occurrence (%)')

# Title
# title = ax.set_title('Average over 2008-2017')

# Add platform locations
# Make sure the order of the platforms corresponds to other notebooks: 
# ['MMIJ', 'LEG', 'BWF1', 'BWF2', 'EPL', 'HKZA', 'HKZB', 'HKNA', 'HKNB', 'K13']
platforms = [['MMIJ', 3.435566666666667, 52.84816666666667, 0, 5],
             ['LEG', 3.6670277777777778, 51.91711111111111, 2, -.5],
             ['BWF1', 3.034638888888889, 51.70688888888889, 2, -1],
             ['BWF2', 2.951416666666667, 51.64630555555556, -3.5, -.7],
             ['EPL', 3.276, 51.999, -3.5, 0],
             ['HKZA', 4.008972222222222, 52.30658333333333, 3, .5],
             ['HKZB', 4.008583333333333, 52.28908333333333, 3, 0],
             ['HKNA', 4.242020474993242, 52.68869454787769, 1, 3],
             ['HKNB', 4.241912343238767, 52.683319381733675, 1.8, 2.2],
             ['K13', 3.218611111111111, 53.2175, -2, 3]]

transform = ccrs.PlateCarree()._as_mpl_transform(ax)
for name, lon, lat, offset_x, offset_y in platforms:
    scat = ax.scatter(lon, lat, transform=ccrs.PlateCarree())
    col = scat.get_facecolors()[0].tolist()
    ax.annotate(name, (lon, lat), (lon+offset_x, lat+offset_y), xycoords=transform, fontsize=11, color=col,
                arrowprops=dict(arrowstyle='->', color=col, linewidth=2.5))
    
ax0.set_title('A. Conceptual illustration of LLJ')
ax.set_title('B. Preliminary spatial LLJ distribution')

fig.savefig('../../reports/figures/llj_illustration.png', bbox_inches='tight')
plt.show()

LLJ maps grouped by lamb weather type

In [22]:
percentages = lwt.value_counts()/len(lwt)*100
percentages.head()
Out[22]:
A     20.689616
W     11.092481
SW     9.878867
C      8.275162
NW     6.110275
Name: LWT, dtype: float64
In [27]:
# Now groupby lamb weather type
lljs_by_lwt = lljs.groupby(lwt.to_xarray().rename({'index':'time'})).mean(dim='time').sel(LWT=mapping)*100
msl_by_lwt = msl.msl.groupby(lwt.to_xarray().rename({'index':'time'})).mean(dim='time').sel(LWT=mapping)/100 # to hPa
lon = lljs.longitude
lat = lljs.latitude
lon_msl = msl.longitude
lat_msl = msl.latitude

# Use a lambert projection centered around the IJmuiden mast
clat,clon = 52.8481667, 3.4356667
map_proj = ccrs.LambertConformal(
    central_longitude=clon,
    central_latitude=clat,
    standard_parallels=(45,55),
    globe=ccrs.Globe(datum='WGS84',ellipse='sphere'))

# Set up the figure, use fixed number of clusters for convenience
fig, axs = plt.subplots(9, 3, figsize=(12, 30), dpi=100, subplot_kw={
    'projection': map_proj, 'extent': [-8,15,48,60]})

axs = axs.flatten()
borders = NEF(category='cultural', name='admin_0_countries', scale='50m', facecolor='none') # 110, 50 or 10
for wt, ax in zip(mapping, axs):
    meash = ax.pcolormesh(lon, lat, lljs_by_lwt.sel(LWT=wt).values, vmin=0, vmax=30, cmap='GnBu', transform=ccrs.PlateCarree())
    ax.outline_patch.set_visible(False)
    ax.set_title(wt+'  ({:.1f}%)'.format(percentages[wt]), fontsize=18)
    ax.add_feature(borders, edgecolor='gray')
    
    u, v = geowind(msl_by_lwt.sel(LWT=wt))
    speed = (u**2+v**2)**.5
    lw = 2*speed/speed.max()
    ax.streamplot(lon_msl, lat_msl, u, v, color='grey', linewidth=lw, transform=ccrs.PlateCarree(), density=0.6)

# fig.tight_layout()
plt.colorbar(mesh, ax=axs, label='Average LLJ rate based on 10 years of ERA5 data up to 500m (%)' )

# fig.savefig('../../reports/figures/era5_llj_by_lwt', transparent=False, bbox_inches='tight')
Out[27]:
<matplotlib.colorbar.Colorbar at 0x2aecab7d5630>
In [104]:
# Same horizontal
# LWTs are automatically sorted on groupby, I want to sort them in a different order:
mapping = ["A", "AN", "ANE", "AE", "ASE", "AS", "ASW", "AW", "ANW",
           "C", "CN", "CNE", "CE", "CSE", "CS", "CSW", "CW", "CNW",
           "U",  "N",  "NE",  "E",  "SE",  "S",  "SW",  "W",  "NW"]
In [28]:
# Now groupby lamb weather type
lljs_by_lwt = lljs.groupby(lwt.to_xarray().rename({'index':'time'})).mean(dim='time').sel(LWT=mapping)*100
msl_by_lwt = msl.msl.groupby(lwt.to_xarray().rename({'index':'time'})).mean(dim='time').sel(LWT=mapping)/100 # to hPa
lon = lljs.longitude
lat = lljs.latitude
lon_msl = msl.longitude
lat_msl = msl.latitude

# Use a lambert projection centered around the IJmuiden mast
clat,clon = 52.8481667, 3.4356667
map_proj = ccrs.LambertConformal(
    central_longitude=clon,
    central_latitude=clat,
    standard_parallels=(45,55),
    globe=ccrs.Globe(datum='WGS84',ellipse='sphere'))

# Set up the figure, use fixed number of clusters for convenience
fig, axs = plt.subplots(3, 9, figsize=(30, 11), dpi=100, subplot_kw={
    'projection': map_proj, 'extent': [-8,15,47,63]})

axs = axs.flatten()
borders = NEF(category='cultural', name='admin_0_countries', scale='50m', facecolor='none') # 110, 50 or 10
for wt, ax in zip(mapping, axs):
    mesh = ax.pcolormesh(lon, lat, lljs_by_lwt.sel(LWT=wt).values, vmin=0, vmax=25, cmap='GnBu', transform=ccrs.PlateCarree())
    ax.outline_patch.set_visible(False)
    ax.set_title(wt+'  ({:.1f}%)'.format(percentages[wt]), fontsize=18)
    ax.add_feature(borders, edgecolor='gray')
    
    u, v = geowind(msl_by_lwt.sel(LWT=wt))
    speed = (u**2+v**2)**.5
    lw = 2*speed/speed.max()
    ax.streamplot(lon_msl, lat_msl, u, v, color='grey', linewidth=lw, transform=ccrs.PlateCarree(), density=0.6)

cbar_ax = fig.add_axes([0.15, 0.05, 0.7, 0.02])
cbar = plt.colorbar(mesh, cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=18)
cbar.set_label('Average LLJ rate based on 10 years of ERA5 data up to 500m (%)', fontsize=18)

fig.savefig('../../reports/figures/era5_llj_by_lwt_hor', transparent=False, bbox_inches='tight')

Let's try some different layouts that might fit better ???

In [29]:
clat,clon = 52.8481667, 3.4356667
map_proj = ccrs.LambertConformal(
    central_longitude=clon,
    central_latitude=clat,
    standard_parallels=(45,55),
    globe=ccrs.Globe(datum='WGS84',ellipse='sphere'))

# Set up the figure, use fixed number of clusters for convenience
fig = plt.figure(figsize=(30, 20), dpi=100)
axs = []
for relpos in range(5):
    axs.append(fig.add_axes([relpos/5+1/10, 5/6, 1/5.1, 1/6.1]))
    axs.append(fig.add_axes([relpos/5, 4/6, 1/5.1, 1/6.1]))
    axs.append(fig.add_axes([relpos/5+1/10, 3/6, 1/5.1, 1/6.1]))
    axs.append(fig.add_axes([relpos/5, 2/6, 1/5.1, 1/6.1]))
    axs.append(fig.add_axes([relpos/5+1/10, 1/6, 1/5.1, 1/6.1]))
    axs.append(fig.add_axes([relpos/5, 0/6, 1/5.1, 1/6.1]))
    
plt.show()
In [30]:
fig, axs = plt.subplots(3,3)
In [31]:
fig = plt.figure()
k=0
for i in range(3):
    for j in range(8):
        k+=1
        fig.add_subplot(4,8,k)
        
fig.add_subplot(4,8,26)
fig.add_subplot(4,8,28)
fig.add_subplot(4,8,30)
Out[31]:
<matplotlib.axes._subplots.AxesSubplot at 0x2aecb1e71400>
In [32]:
fig = plt.figure()
k=8
for i in range(3):
    for j in range(8):
        k+=1
        fig.add_subplot(4,8,k)
        
fig.add_subplot(4,8,2)
fig.add_subplot(4,8,4)
fig.add_subplot(4,8,6)
Out[32]:
<matplotlib.axes._subplots.AxesSubplot at 0x2aecb246e198>
In [33]:
# monkey patch to fix incompatibility between cartopy and matplotlib 
# https://github.com/SciTools/cartopy/issues/1120#issuecomment-424418760
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh

# Same horizontal
# LWTs are automatically sorted on groupby, I want to sort them in a different order:
mapping = ["A", "U", "C",
           "AN", "ANE", "AE", "ASE", "AS", "ASW", "AW", "ANW",
           "N",  "NE",  "E",  "SE",  "S",  "SW",  "W",  "NW",
           "CN", "CNE", "CE", "CSE", "CS", "CSW", "CW", "CNW"]

# Use a lambert projection centered around the IJmuiden mast
clat,clon = 52.8481667, 3.4356667
map_proj = ccrs.LambertConformal(
    central_longitude=clon,
    central_latitude=clat,
    standard_parallels=(45,55),
    globe=ccrs.Globe(datum='WGS84',ellipse='sphere'))

# Set up the figure
fig = plt.figure(figsize=(30, 16), dpi=100)
axs = []
axs.append(fig.add_subplot(4,8,2, projection=map_proj, extent=[-8,15,47,63]))
axs.append(fig.add_subplot(4,8,4, projection=map_proj, extent=[-8,15,47,63]))
axs.append(fig.add_subplot(4,8,6, projection=map_proj, extent=[-8,15,47,63]))

k=8
for i in range(3):
    for j in range(8):
        k+=1
        axs.append(fig.add_subplot(4,8,k, projection=map_proj, extent=[-8,15,47,63]))       

# Now groupby lamb weather type
lljs_by_lwt = lljs.groupby(lwt.to_xarray().rename({'index':'time'})).mean(dim='time').sel(LWT=mapping)*100
msl_by_lwt = msl.msl.groupby(lwt.to_xarray().rename({'index':'time'})).mean(dim='time').sel(LWT=mapping)/100 # to hPa
lon = lljs.longitude
lat = lljs.latitude
lon_msl = msl.longitude
lat_msl = msl.latitude

# Plot the data
borders = NEF(category='cultural', name='admin_0_countries', scale='50m', facecolor='none') # 110, 50 or 10
for wt, ax in zip(mapping, axs):
    mesh = ax.pcolormesh(lon, lat, lljs_by_lwt.sel(LWT=wt).values, vmin=0, vmax=25, cmap='GnBu', transform=ccrs.PlateCarree())
    ax.outline_patch.set_visible(False)
    ax.set_title(wt+'  ({:.1f}%)'.format(percentages[wt]), fontsize=18)
    ax.add_feature(borders, edgecolor='gray')
    
    u, v = geowind(msl_by_lwt.sel(LWT=wt))
    speed = (u**2+v**2)**.5
    lw = 2*speed/speed.max()
    ax.streamplot(lon_msl, lat_msl, u, v, color='grey', linewidth=lw, transform=ccrs.PlateCarree(), density=0.6)

cbar_ax = fig.add_axes([0.15, 0.05, 0.7, 0.02])
cbar = plt.colorbar(mesh, cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=18)
cbar.set_label('Average LLJ rate based on 10 years of ERA5 data up to 500m (%)', fontsize=18)

fig.savefig('../../reports/figures/era5_llj_by_lwt_alt', transparent=False, bbox_inches='tight')
fig.savefig('../../reports/figures/era5_llj_by_lwt_alt_dpi200', transparent=False, bbox_inches='tight', dpi=200)