In [2]:
location = 'Struisbaai'
In [1]:
import numpy as np
import xarray as xr
import pandas as pd
import glob

from ecmwf.opendata import Client
client = Client()

from datetime import datetime, timedelta

import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.express as px

import math
import pickle
import plotly.io as pio

Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.16.0
In [3]:
today = datetime.today().date()
In [4]:
for i in range(3):
    date = str(today - timedelta(days = i))
    print(date)
    client.retrieve(
     date=str(date),
     time=0,
     step=list(np.arange(0,24,3)),
     stream="enfo",
     type=['cf'],
     levtype="sfc",
     param=["10u","10v"],
     target='ECMWF_forecast/10m_U_V_'+str(date)+'.grib'
    )

2022-04-13

<multiple>:   0%|          | 0.00/9.29M [00:00<?, ?B/s]
2022-04-12

<multiple>:   0%|          | 0.00/9.29M [00:00<?, ?B/s]
2022-04-11

<multiple>:   0%|          | 0.00/9.29M [00:00<?, ?B/s]
In [5]:
client.retrieve(
 date=str(today),
 time=0,
 step=list(np.arange(24,147,3)),
 stream="enfo",
 type=['cf'],
 levtype="sfc",
 param=["10u","10v"],
 target='ECMWF_forecast/10m_U_V_7day_forecast.grib'
)

<multiple>:   0%|          | 0.00/47.6M [00:00<?, ?B/s]
Out[5]:
<ecmwf.opendata.client.Result at 0x7f59f41c4f70>
In [4]:
files = glob.glob('ECMWF_forecast/10m_U_V_*.grib')
files.sort()
In [5]:
grbs = []
for file in files:
    ds = xr.open_dataset(file)
    ds = ds.drop('time')
    ds = ds.rename_dims({"step":"time"}) 
    ds['time'] = ds.valid_time
    grbs.append(ds)
    
ds = xr.merge(grbs)

Ignoring index file 'ECMWF_forecast/10m_U_V_2022-04-11.grib.923a8.idx' older than GRIB file
Ignoring index file 'ECMWF_forecast/10m_U_V_2022-04-12.grib.923a8.idx' older than GRIB file
Ignoring index file 'ECMWF_forecast/10m_U_V_7day_forecast.grib.923a8.idx' older than GRIB file
In [6]:
ds['u10'] = ds.u10*1.94384
ds['v10'] = ds.v10*1.94384
In [7]:
def wind_speed(U, V):
    func = lambda x, y: np.sqrt(x**2 + y**2)
    return xr.apply_ufunc(func, U, V, dask = 'parallelized')

def wind_dir(U, V):
    func = lambda x, y: np.mod(180+np.rad2deg(np.arctan2(x, y)),360)
    return xr.apply_ufunc(func, U, V, dask = 'parallelized')
In [8]:
ds['ws'] = wind_speed(ds.u10, ds.v10)
ds['dir'] = wind_dir(ds.u10, ds.v10)
In [9]:
latlong_dict = {}
latlong_dict['Struisbaai'] = [-34.80, 20.10]
latlong_dict['Mossel Bay'] = [-34.15, 22.20]
latlong_dict['Muizenberg'] = [-34.11, 18.52]
latlong_dict['Witsands'] = [-34.41, 20.92]
latlong_dict['Arniston'] = [-34.67, 20.27]
latlong_dict['Plettenberg Bay'] = [-34.07, 23.44]
latlong_dict['Hermanus'] = [-34.44, 19.25]
latlong_dict['Sedgefield'] = [-34.10, 22.78]
latlong_dict['Gaansbaai'] = [-34.64, 19.35]

coords = latlong_dict[location]
In [10]:
ds = ds.sel(latitude = coords[0], method = 'nearest').sel(longitude = coords[0],method = 'nearest')
In [13]:
ds_daily = ds.resample(time = 'D', closed = "right").mean().load()

df = ds_daily.to_dataframe()

df['index'] = pd.to_datetime(df.index)

u10_1 = list(ds_daily.u10.values[1:])
u10_1.append(np.nan)
df['u10_1'] = u10_1 
df['u10_5m'] = ds_daily.u10.rolling(time = 5).mean().values
df['month'] = pd.to_datetime(ds_daily.time.values).month

df_ = df[['month', 'ws', 'u10_1', 'u10_5m']].dropna()

features = np.array(df_)

from sklearn.ensemble import RandomForestRegressor

rf = pickle.load(open('model_no_sst_full.pkl', 'rb'))
df_['predict'] = rf.predict(features)

df['predict'] = df_.predict

df['predict'] = df['predict'].fillna(-0.5) # set nan to -0.5 to show days with insufficient data
In [14]:
def get_vector(direction):
    if direction  < 90 and direction > 0:
        x = math.sin(math.radians(direction))
        y = math.cos(math.radians(direction))
    elif direction > 90 and direction <180:
        x = math.sin(math.radians(180 - direction))
        y = -math.cos(math.radians(180 - direction))
    elif direction >180 and direction <270:
        x = -math.cos(math.radians(270 - direction))
        y = -math.sin(math.radians(270 - direction))
    else:
        x = -math.sin(math.radians(360-direction))
        y = math.cos(math.radians(360-direction))
    return(x,-y)
In [15]:
times = pd.to_datetime(ds.time.values)
arrows = []
for dir_,i in zip(ds.dir.values[0::2], times[0::2]):
    ax_,ay_ = get_vector(dir_)
    arrows.append(dict(
        x=i,
        y=-5,
        xref = 'x',
        yref = 'y',
        axref = 'pixel',
        ayref = 'pixel',
        ax=ax_*25,
        ay=ay_*25,
        arrowhead = 2,
        arrowsize = 1,
        arrowwidth = 1.4,
        xanchor = 'left',
        yanchor = 'bottom',
        visible = True
    )
    )
In [16]:
fig = make_subplots(rows=2, cols=1, row_heights=[0.7, 0.3])

fig.add_trace(
    go.Scatter(x = times, y= ds.ws, name='Wind Speed (knots)', visible=True), row = 1, col =1)


x = df.predict.where(df.predict<0).dropna()
fig.add_trace(
    go.Scatter(x = x.index.values, y= list(x.values),  
                    visible=True, 
                    marker = dict(size=10, color = 'grey'),
                    mode='markers',
                    name='Insufficent data'), row = 2, col =1)


x = df.predict.where(df.predict>=0).where(df.predict<0.6).dropna()
fig.add_trace(
    go.Scatter(x = x.index.values, y= list(x.values),  
                    visible=True, 
                    marker = dict(size=10, color = 'green'),
                    mode='markers',
                    name='Stranding Unlikely'), row = 2, col =1)

x = df.predict.where(df.predict>=0.6).where(df.predict<1).dropna()
fig.add_trace(
    go.Scatter(x = x.index.values, y= list(x.values),  
                    visible=True, 
                    marker = dict(size=10, color = 'yellow'),
                    mode='markers',
                    name='Stranding Possible'), row = 2, col =1)


x = df.predict.where(df.predict>=1).dropna()
fig.add_trace(
    go.Scatter(x = x.index.values, y= list(x.values),  
                    visible=True, 
                    marker = dict(size=10, color = 'red'),
                    mode='markers',
                    name='Stranding Likely'), row = 2, col =1)



for a in arrows:
    fig.add_annotation(a)
    
fig.update_yaxes(range=[-10, 30], row = 1)
fig.update_xaxes(range=[times[-64], times[-1]], row = 1)

if df.predict.max() < 1:
    fig.update_yaxes(range=[-1, 5], row = 2)
    
fig.update_xaxes(range=[times[-64], times[-1]], row = 2)


fig.update_yaxes(title="Wind Speed(knots)", row = 1)
fig.update_yaxes(title="Strandings Projected", row = 2)

fig.show()
In [17]:
pio.write_html(fig, file='forecast.html', auto_open=False)
In [ ]:
In [ ]:
In [ ]:

end