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 (1/6), we provide a brief introduction into the observation data.
The observation data come from 7 sites. LiDARs were installed at each site, and three sites were equipped with two LiDARs. These are the abbreviations.
A website dedicated to these measurement is available at https://windopzee.net/en/home/. The data are stored as NetCDF format, 1 file per platform:
ls ../../data/external/ECN/*.nc
The postprocessing of the platform data is described in the appendix of the paper. The resulting output files all have the following format.
import xarray as xr
ds = xr.open_dataset('../../data/external/ECN/BWFZwinds_lot2.nc')
ds
The output files were produced using matlab, while the further analysis in done in python. To achieve good compatibility, I wrote a preprocessing function, which addresses some remaining inconsistencies.
def clean(ds):
""" This function modifies the platform datasets on load to addresses the following:
- Convert the units of the time arrays in order to read them properly
- Change variable names (remove spaces etc.)
- Discard UTC time and model reference time (not necessary).
- Set time and height variables as coordinates rather than variables
- Re-index time array to have uniform spacing (convert gaps to NaN values)
"""
import xarray as xr
from pandas import date_range
# Convert units of time arrays (convert capitals to lower case)
ds['Model Reference Time'].attrs['units'] = ds['Model Reference Time'].attrs['units'].lower()
ds['File Reference Time'].attrs['units'] = ds['File Reference Time'].attrs['units'].lower()
# Rename variables
ds.rename({
"File Reference Time": "time",
"Measurement Levels": "height",
"Site Latitude": "latitude",
"Site Longitude": "longitude",
"WS": "wspd",
"WD": "wdir"}, inplace=True)
# Discard redundant time arrays
ds = ds.drop(['Model Reference Time', 'UTC Time'])
# This last function correctly parses the time arrays and makes the time equidistant (filling up with nan)
ds = xr.decode_cf(ds) #.resample(time='10min').asfreq() # Equivalent function, but slower (I think)
return ds.reindex(time=date_range(ds.time.values[0], ds.time.values[-1], freq='10min'))
ds = xr.open_dataset('../../data/external/ECN/BWFZwinds_lot1.nc')
clean(ds)
We now load all the datasets into one dictionary. This code is copied to the code base for later re-use (code/datasets/ecn.py)
data_path = '../../data/external/ECN/'
files = {
'BWF1': data_path+'BWFZwinds_lot1.nc',
'BWF2': data_path+'BWFZwinds_lot2.nc',
'EPL': data_path+'EPLwinds.nc',
'HKNA': data_path+'HKNAwinds.nc',
'HKNB': data_path+'HKNBwinds.nc',
'HKZA': data_path+'HKZAwinds.nc',
'HKZB': data_path+'HKZBwinds.nc',
'K13': data_path+'K13winds.nc',
'LEG': data_path+'LEGwinds.nc',
'MMIJ': data_path+'MMIJwinds.nc'}
# Save the platform data in a dictionary
obs = {}
for name, file in files.items():
print(name)
obs[name] = clean(xr.open_dataset(file))
import matplotlib.pyplot as plt
plt.xkcd()
fig, ax = plt.subplots()
for name, ds in obs.items():
ax.plot(ds.wspd.mean(dim='time').values, ds.height.values, 'o-', label=name)
ax.legend(bbox_to_anchor=(1.05, 1))
ax.set_xlabel('mean wind speed (m/s)')
ax.set_ylabel('altitude (m)')
ax.set_title('Vertical wind profiles for each of the LiDARs')
plt.show()
fig, ax = plt.subplots()
for name, ds in obs.items():
ax.plot(range(24), ds.wspd.mean(dim='height').groupby('time.hour').mean(dim='time'), label=name)
ax.legend(bbox_to_anchor=(1.05, 1))
ax.set_ylabel('height-/time-averaged wind speed (m/s)')
ax.set_xlabel('hour of the day')
ax.set_title('diurnal cycle of wind speed')
plt.show()
The figure below shows the time/height evolution of wind speed for each dataset. It illustrates the temporal extent and overlap of the datasets, as well as the missing data.
fig, axs = plt.subplots(11, 1, figsize=(12, 16), sharex=True, sharey=True)
axs = axs.flatten()
for i, (name, ds) in enumerate(obs.items()):
mesh = axs[i].pcolormesh(ds.time, ds.height, ds.wspd, vmin=0, vmax=30)
axs[i].set_title(name)
plt.tight_layout()
axs[-2].xaxis.set_tick_params(labelbottom=True)
cax = fig.add_axes(axs[-1].get_position())
axs[-1].set_axis_off()
plt.colorbar(mesh, cax=cax, orientation='horizontal', label='wind speed (m/s)')
plt.show()
For wind direction, I want a cyclic colormap The code below is not relevant for understanding the data, but is necessary for the subsequent visualization.
# For wind direction, I want a cyclic colormap
# Use function from great answer: https://stackoverflow.com/a/34557535/6012085
import numpy as np
import matplotlib.colors as col
!pip install hsluv
import hsluv
##### generate custom colormaps
def make_segmented_cmap():
white = '#ffffff'
black = '#000000'
red = '#ff0000'
blue = '#0000ff'
anglemap = col.LinearSegmentedColormap.from_list(
'anglemap', [black, red, white, blue, black], N=256, gamma=1)
return anglemap
def make_anglemap( N = 256, use_hpl = True ):
h = np.ones(N) # hue
h[:N//2] = 11.6 # red
h[N//2:] = 258.6 # blue
s = 100 # saturation
l = np.linspace(0, 100, N//2) # luminosity
l = np.hstack( (l,l[::-1] ) )
colorlist = np.zeros((N,3))
for ii in range(N):
if use_hpl:
colorlist[ii,:] = hsluv.hpluv_to_rgb( (h[ii], s, l[ii]) )
else:
colorlist[ii,:] = hsluv.hsluv_to_rgb( (h[ii], s, l[ii]) )
colorlist[colorlist > 1] = 1 # correct numeric errors
colorlist[colorlist < 0] = 0
return col.ListedColormap( colorlist )
N = 256
hpluv_anglemap = make_anglemap( use_hpl = True )
The plot below is similar to the previous one for wind speed, but now for wind direction. This is not used in the paper so it is less relevant. In this case, the time axes are not shared between the axis, so it is in fact a 'zoomed in' version of the previous plot.
fig, axs = plt.subplots(11, 1, figsize=(12, 16), sharex=False, sharey=False)
axs = axs.flatten()
for i, (name, ds) in enumerate(obs.items()):
mesh = axs[i].pcolormesh(ds.time, ds.height, ds.wdir, vmin=0, vmax=360, cmap=hpluv_anglemap)
axs[i].set_title(name)
plt.tight_layout()
axs[-2].xaxis.set_tick_params(labelbottom=True)
cax = fig.add_axes(axs[-1].get_position())
axs[-1].set_axis_off()
plt.colorbar(mesh, cax=cax, orientation='horizontal', label='wind direction')
plt.show()
for name, ds in obs.items():
ds.to_netcdf('../../data/interim/ECNplatforms/'+name+'.nc')
So far so good! Some of the scripts used in this notebook have been refactored to a shared code base in order to load the data more quickly in other notebooks.