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:
In this notebook (6/6), I create some of the spatial plots that are based on large datasets of ERA5 data.
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()
# 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']
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
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
# 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()
percentages = lwt.value_counts()/len(lwt)*100
percentages.head()
# 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')
# 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"]
# 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')
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()
fig, axs = plt.subplots(3,3)
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)
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)
# 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)