Contents

MÉRA-based Wind Atlas for Irish Continental Shelf region

_images/eooffshore_banner.png 

_images/seai.png _images/ucd.png

MÉRA-based Wind Atlas for Irish Continental Shelf region

Wind Atlas Prototype

Wind atlases have been developed to provide energy resources maps, containing information on wind speeds and related variables at multiple heights above sea level for offshore areas of interest (AOIs). One example that focuses on Irish offshore AOIs is the Sustainable Energy Authority of Ireland (SEAI) Wind Atlas, which is described by UK Met Office, 2015 - Remodelling the Irish national onshore and offshore wind atlas.

This notebook demonstrates the utility of the Pangeo software ecosystem in the development of an Irish offshore wind atlas prototype, covering offshore renewable energy assessment areas in the Irish Continental Shelf (ICS) region. It uses the EOOffshore catalog MÉRA data set created for this region, which is an analysis-ready, cloud-optimized (ARCO) dataset featuring 16 years of reanalysis wind data products. Scalable processing and visualization of this ARCO catalog is demonstrated with analysis of provided data variables and computation of new variables as required for AOIs, avoiding redundant storage and processing requirements for areas not under assessment. The MÉRA data set is described in the MÉRA Wind Data for Irish Continental Shelf region notebook.

The prototype implementation provides:

How to cite: O’Callaghan, D. and McBreen, S.: Scalable Offshore Wind Analysis With Pangeo, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-2746, https://doi.org/10.5194/egusphere-egu22-2746, 2022.

Note:

  • Full interactive dashboard functionality is only available when this notebook is executed.

  • If viewing this notebook on the EOOffshore website:

    • Click the “Fullscreen Mode” button at the top of the page to ensure that the wind atlas is correctly displayed.

    • The dashboard will be displayed with limited functionality.

    • Code cells have been hidden, but may be viewed via the corresponding “> Click to show” buttons.

%matplotlib inline
from dataclasses import dataclass, field
from datetime import datetime
import geopandas as gpd
from intake import Catalog, open_catalog
import json
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import regionmask
import seaborn as sns
import shapely.geometry as sgeom
from sklearn.neighbors import NearestNeighbors
from typing import Dict, List, Tuple
from windrose import WindroseAxes
import xarray as xr

import holoviews as hv
hv.extension("bokeh")
import datashader as dsh
import geoviews as gv
import geoviews.feature as gf
from holoviews.operation.datashader import rasterize
import hvplot.pandas
import hvplot.xarray
import panel as pn

import warnings
warnings.filterwarnings('ignore')

catalog = open_catalog('data/intake-catalogs/eooffshore_ics.yaml')

# Wind turbine speed cut-in/cut-out limits (m/s)
CUT_IN_WIND_SPEED = 3
CUT_OUT_WIND_SPEED = 25

@dataclass
class EooffshoreDataset:
    """Base EOOffshore dataset class, associated with a single Intake Catalog entry"""
    name: str # Short name used in plots
    catalog: Catalog = field(repr=False, hash=False, compare=False)
    dataset_id: str # Catalog identifier

@dataclass
class XarrayEooffshoreDataset(EooffshoreDataset):
    """Dataset backed by a Zarr store, NetCDF file etc."""
    dataset: xr.Dataset = field(init=False, hash=False, compare=False, repr=False)
    rename_variables: dict[str, str] = None
    rotate_longitude: bool = False # Rotate from [0, 360] to [-180, 180]
    
    def __post_init__(self):
        """Load the dataset from the Intake catalog, and perform some common initialisation"""
        self.dataset = self.catalog[self.dataset_id].to_dask()
        
        # Enable consistent variable naming
        if self.rename_variables:
            self.dataset = self.dataset.rename(self.rename_variables)
           
    def sel_grid_nearest(self, latitude, longitude):
        """Select a subset of the dataset, nearest to the specified coordinates"""
        if self.rotate_longitude:
            longitude %= 360     
        return self._nearest_latitude_longitude(latitude, longitude)
            
    def _nearest_latitude_longitude(self, latitude, longitude):
        """Select a subset of the associated 1-D grid dataset, nearest to the specified coordinates"""
        return self.dataset.sel(longitude=longitude, latitude=latitude, method='nearest')

@dataclass
class CurvilinearEooffshoreDataset(XarrayEooffshoreDataset):
    """
    Additional functionality to handle issue where `xarray.Dataset.sel()` doesn't support calls to sel() for 2-D latitude/longitude grids. 
    This is enabled by the use of a K-Nearest Neighbours (kNN) learner
    """
    grid_neighbours: NearestNeighbors = field(init=False, hash=False, compare=False, repr=False)

    def __post_init__(self):
        """Fit nearest neighbours learner using dataset latitude/longitude grid, used in subsequent calls to sel()"""
        super().__post_init__()
        
        longitude = self.dataset.longitude
        latitude = self.dataset.latitude

        self.grid_neighbours = NearestNeighbors(n_neighbors=1, algorithm='auto').fit(np.stack((longitude.values.reshape(-1,), latitude.values.reshape(-1,)), axis=1))

    def _nearest_latitude_longitude(self, latitude, longitude):
        """Select a subset of the associated 2-D grid dataset, using the top kNN (k=1) neighbours to the specified coordinates"""
        # Select one neighbour similar to `xarray.Dataset.sel(..., method='nearest')`
        k = 1
        distances, indices = self.grid_neighbours.kneighbors(np.array([longitude, latitude]).reshape(1,-1), n_neighbors=k)
        yindex, xindex = np.unravel_index(indices, self.dataset.latitude.shape)
        ydim, xdim = self.dataset.latitude.dims
        return self.dataset.isel({ydim: yindex[0][0], xdim: xindex[0][0]})
    
@dataclass
class WindfarmEooffshoreDataset(EooffshoreDataset):
    """Wind farm location shapefile dataset"""
    dataframe: gpd.GeoDataFrame = field(init=False, hash=False, compare=False, repr=False)
    windfarms: List[str] = None 
    windfarm_mask: xr.DataArray = field(init=False, hash=False, compare=False, repr=False)

    def __post_init__(self):
        """Load the wind farm dataset from the catalog shapefile, and create a corresponding mask used for subsequent dataset subset selection"""
        self.dataframe = gpd.read_file(self.catalog[self.dataset_id].urlpath).to_crs("EPSG:4326")
        self.windfarms = self.dataframe.Name.to_list()
        
        # Higher-resolution mask that's used when dataset mask is insufficient during subsequent windfarm AOI retrieval
        self.windfarm_mask = regionmask.mask_geopandas(self.dataframe, np.arange(-11, -5, .007), np.arange(55.42, 51.0, -.007)).rename({'lat': 'latitude', 'lon':'longitude'})
    
@dataclass
class Buoy:
    name: str
    longitude: float
    latitude: float
    
@dataclass
class IwbEooffshoreDataset(EooffshoreDataset):
    """Irish Weather Buoy dataset"""
    dataframe: pd.DataFrame = field(init=False, hash=False, compare=False, repr=False)
    buoy_coordinates: Dict[str, Buoy] = field(init=False)
    
    def __post_init__(self):
        """Preprocessing of the data frame loaded from the Intake catalog"""
        self.dataframe = self.catalog[self.dataset_id].read()
        # Retain only buoys still in use
        self.dataframe = self.dataframe[self.dataframe.station_id.isin(['M2', 'M3', 'M4', 'M5', 'M6'])].rename(columns={'longitude (degrees_east)': 'longitude', 'latitude (degrees_north)': 'latitude'})
        self.buoy_coordinates = {name: Buoy(name=name, longitude=x['longitude'], latitude=x['latitude']) for name, x in self.dataframe.groupby('station_id')[['longitude', 'latitude']].agg(['unique']).applymap(lambda x:x[0]).droplevel(1, axis=1).T.to_dict().items()}

# ERA5 data used for map plot limits
era5 = XarrayEooffshoreDataset(name='ERA5', catalog=catalog, dataset_id='eooffshore_ics_era5_single_level_hourly_wind')

mera = CurvilinearEooffshoreDataset(name='MÉRA', 
                                    catalog=catalog, 
                                    dataset_id='eooffshore_ics_mera_wind_ANALYSIS', 
                                    rotate_longitude=True,
                                   )
mera.dataset.atlas_mean_wind_speed.attrs = mera.dataset.wind_speed.attrs

# Some variables have been precomputed for use with interactive dashboard
WIND_ATLAS_VARIABLES = ['atlas_mean_wind_speed', 
                        'atlas_mean_power_density',
                        'atlas_weibull_shape',
                        'atlas_weibull_scale',
                        'atlas_operational_frequency',
                        'atlas_maximum_yield_frequency']
for variable in WIND_ATLAS_VARIABLES + ['latitude', 'longitude']:
    mera.dataset[variable] = mera.dataset[variable].compute()

# IWB buoys
iwb = IwbEooffshoreDataset(name='IWB', catalog=catalog, dataset_id='eooffshore_iwb_hourly')

# Synthetic windfarms
eoowindfarmds = WindfarmEooffshoreDataset(name='EOOffshore WF', catalog=catalog, dataset_id='eooffshore_windfarms')

# Load the Irish Continental Shelf geometry
with open('data/linestring_ics_geo.json') as jf:
    icsgeo = json.load(jf)
icsline = sgeom.LineString(icsgeo['features'][0]['geometry']['coordinates'])

WINDFARM_COLORS = [sns.color_palette('tab20')[x] for x in [9, 3, 10, 5, 8]]

lim_buffer = 0.05
ICS_XLIM = (era5.dataset.longitude.min().item() - lim_buffer, era5.dataset.longitude.max().item() + lim_buffer)
ICS_YLIM = (era5.dataset.latitude.min().item() - lim_buffer, era5.dataset.latitude.max().item() + lim_buffer)

# The windrose projection doesn't render bars correctly unless at least one polar projection plot has been created
# For now, this placeholder appears to be sufficient for subsequent correct rendering
N = 1
theta = np.linspace(0.0, 2 * np.pi, N, endpoint=False)
radii = 10 * np.random.rand(N)
width = np.pi / 4 * np.random.rand(N)
ax = plt.subplot(111, projection='polar')
bars = ax.bar(theta, radii, width=width, bottom=0.0)
plt.close()
# <end placeholder>

# https://github.com/numpy/numpy/issues/18363#issuecomment-775101347
format_plot_time = lambda time: time.values.astype('datetime64[us]').astype(datetime).strftime('%Y/%m')

# Map longitude must be in [-180, 180] to ensure zoom works
power_map_ds = mera.dataset[['atlas_mean_power_density']].rename({'latitude': 'Latitude', 'longitude': 'Longitude', 'height': 'Height', 'atlas_mean_power_density': 'Power Density'})
power_map_ds['Longitude'] = power_map_ds.Longitude - 360
power_map = rasterize(gv.Dataset(data=power_map_ds, 
                                 kdims=['Height', 'Longitude', 'Latitude'], 
                                 vdims='Power Density').to(gv.QuadMesh, ['Longitude', 'Latitude']), precompute=True).options(clim=(400, 1700), 
                                                                                                                                   cmap=plt.cm.get_cmap('inferno_r', 15), 
                                                                                                                                   width=590, 
                                                                                                                                   height=450, 
                                                                                                                                   colorbar=True, 
                                                                                                                                   tools=['hover'],
                                                                                                                                   title=f'Power Density ({format_plot_time(mera.dataset.time.min())} - {format_plot_time(mera.dataset.time.max())})',
                                                                                                                                   fontsize={'title': 9},
                                                                                                                                   )

mapstream = hv.streams.Tap(source=power_map, x=-12, y=53)

def selected_height():
    selected_height = power_map.dimension_values('Height')
    if len(selected_height) == 0: # Upon creation
        selected_height = mera.dataset.height.values[0]
    return int(selected_height)

@pn.depends(mapstream.param.x, mapstream.param.y)
def windspeed_timeseries(x, y):
    nearest = mera.sel_grid_nearest(latitude=y, longitude=x)
    height = selected_height()
    ws = nearest.sel(height=height).wind_speed.resample(time='M').mean().rename('Wind speed').rename({'time': 'Time'})
    return ws.hvplot('Time', 
                     title=f'{height}m Wind Speed (Monthly Mean)',
                     fontsize={'title': 9},
                     ylim=(CUT_IN_WIND_SPEED, 19), 
                     ylabel=f'{nearest.wind_speed.long_name} [{nearest.wind_speed.units}]'
                    )

@pn.depends(mapstream.param.x, mapstream.param.y)
def windspeed_month(x, y):
    nearest = mera.sel_grid_nearest(latitude=y, longitude=x)
    height = selected_height()
    ws = nearest.sel(height=height).wind_speed.groupby('time.month').mean().rename('Wind speed').rename({'month': 'Month'})
    return ws.hvplot('Month', 
                       title=f'Mean {height}m Wind Speed (Month of Year)',
                       fontsize={'title': 9},
                       ylim=(CUT_IN_WIND_SPEED, 19), 
                       width=350, 
                       ylabel=f'{nearest.wind_speed.long_name} [{nearest.wind_speed.units}]'
                      ).opts(default_tools=[])

@pn.depends(mapstream.param.x, mapstream.param.y)
def windspeed_dist(x, y):
    nearest = mera.sel_grid_nearest(latitude=y, longitude=x)
    return nearest.sel(height=selected_height()).hvplot.hist(y='wind_speed', xlim=(-0.5, 30), width=330).opts(default_tools=[])

@pn.depends(mapstream.param.x, mapstream.param.y)
def coordinate_statistics(x,y):
    # https://github.com/holoviz/holoviews/issues/3918#issuecomment-523865218
    def hide_index(plot, element):
        plot.handles['table'].index_position = None

    height = selected_height()
    nearest = mera.sel_grid_nearest(latitude=y, longitude=x)
    df = nearest.sel(height=slice(height, height))[WIND_ATLAS_VARIABLES].to_dataframe().reset_index().T
    format_ints = ['height', 'atlas_mean_power_density']
    for x in df.index:
        precision = 0 if x in format_ints else 2
        df.loc[x] = f'{df.loc[x].values[0]:,.{precision}f}'
    df.index = [f'{nearest[variable].long_name} [{nearest[variable].units}]'.capitalize() for variable in df.index]

    return df.reset_index().rename(columns={'index': ' ', 0: '(Mean) Value'}).hvplot.table(width=410, 
                                                                                           height=250, 
                                                                                           title='Selected Coordinate',
                                                                                           fontsize={'title': 9},
                                                                                          ).opts(hooks=[hide_index])    

@pn.depends(mapstream.param.x, mapstream.param.y)
def wind_profile(x, y):
    nearest = mera.sel_grid_nearest(latitude=y, longitude=x)
    wind_profile_df = nearest[WIND_ATLAS_VARIABLES].to_dataframe().reset_index().rename(columns={'atlas_mean_wind_speed': 'Wind speed (mean)', 'height': 'Height'})
    return wind_profile_df.hvplot.line(x='Wind speed (mean)', 
                                      y='Height', 
                                      xlim=(5.5, 13.5), 
                                      yticks=mera.dataset.height.values,
                                      xlabel=f'{nearest.wind_speed.long_name} [{nearest.wind_speed.units}]',
                                      ylabel=f'{nearest.height.long_name} [{nearest.height.units}]',
                                      width=300).opts(default_tools=[]) * wind_profile_df.hvplot.scatter(x='Wind speed (mean)', y='Height').opts(default_tools=[])

@pn.depends(mapstream.param.x, mapstream.param.y)
def windrose(x, y):
    plot = plt.figure(figsize=(5, 4))
    ax = plot.add_subplot(111, projection="windrose", rmax=25)

    bins = np.arange(0, 25, 5)
    
    nearest = mera.sel_grid_nearest(latitude=y, longitude=x)
    height = selected_height()
    ds = nearest.sel(height=slice(height, height))
    ax.bar(direction=ds.wind_direction.values.reshape(-1,), 
           var=ds.wind_speed.values.reshape(-1,), 
           nsector=12, 
           bins=bins, 
           normed=True, 
           opening=0.8, 
           edgecolor='white',
           cmap=plt.cm.YlGnBu)
    ax.set_legend(decimal_places=0, title=mera.dataset.wind_speed.units, title_fontsize='small', loc='lower right')
    labels = ax.get_legend().get_texts()
    for i in range(0, len(bins)-1):
        labels[i].set_text(f'{bins[i]} - {bins[i+1]}')
    labels[-1].set_text(f'{bins[-1]}+')                   
    # Following example in https://anaconda.org/jbednar/panel_gapminders/notebook
    plt.close(plot)
    return plot

ics = gv.Path([icsgeo['features'][0]['geometry']['coordinates']]).opts(gv.opts.Path(line_width=2))
windfarms = gv.Polygons(eoowindfarmds.dataframe.rename(columns={'Name': 'Windfarm'})).opts(cmap=LinearSegmentedColormap.from_list('eoowindfarmds', WINDFARM_COLORS, N=len(eoowindfarmds.dataframe)), tools=['hover'])
iwb_buoys = gv.Points([(b.longitude, b.latitude, b.name) for b in iwb.buoy_coordinates.values()], kdims=['Longitude', 'Latitude'], vdims=['Buoy']).opts(gv.opts.Points(size=9, color=sns.color_palette()[1], tools=['hover']))

aoi_power_map = (power_map * gf.coastline(scale='10m') * ics * iwb_buoys * windfarms).redim.range(Longitude=ICS_XLIM, Latitude=ICS_YLIM)

description = '### [EOOffshore](https://eooffshore.github.io) Wind Speed and Power Density for Irish Continental Shelf region, ' + \
              'using [Met Éireann Re-Analysis (MÉRA)](https://www.met.ie/climate/available-data/mera) data\n' + \
              'Example Areas Of Interest (AOI) including selected [Irish Weather Buoy Network](http://www.marine.ie/Home/site-area/data-services/real-time-observations/irish-weather-buoy-network-imos) ' + \
              'and synthetic wind farm coordinates, are highlighted on the map.'
header = pn.Row(pn.pane.PNG('../eooffshore_banner.png', width=270), pn.pane.PNG('../seai.png', width=250), pn.pane.PNG('../ucd.png', width=60), pn.layout.Spacer(width=30), pn.Pane(description, width=500))

statistics_col = pn.Column(pn.Row(coordinate_statistics, windspeed_dist), pn.Row(pn.panel(windrose), wind_profile))
overview_row = pn.Row(pn.panel(aoi_power_map, widget_location='top'), statistics_col)
windspeed_row = pn.Row(windspeed_timeseries, windspeed_month)
eoo_col = pn.Column(header, overview_row, windspeed_row)
eoo_col

previous

Wind Atlas Development