# 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)