Regional wind distribution

Introduction

Beaufort Wind Scale
0	1	2	3	4	5	6	7	8	9	10	11	12
Calm	Light Air	Light Breeze	Gentle Breeze	Moderate Breeze	Fresh Breeze	Strong Breeze	Near Gale	Gale	Strong Gale	Storm	Violent Storm	Hurricane Force
Light Winds						High Winds		Gale-force		Storm-force		Hurricane-force
<1 mph <1 knot <0.3 m/s	1–3 mph 1–3 knots 0.3–1.5 m/s	4–7 mph 4–6 knots 1.6–3.3 m/s	8–12 mph 7–10 knots 3.4–5.5 m/s	13–18 mph 11–16 knots 5.5–7.9 m/s	18–24 mph 17–21 knots 8.0–10.7 m/s	25–31 mph 22–27 knots 10.8–13.8 m/s	31–38 mph 28–33 knots 13.9–17.1 m/s	39–46 mph 34–40 knots 17.2–20.7 m/s	47-54 mph 41–47 knots 20.8–24.4 m/s	55–63 mph 48–55 knots 24.5–28.4 m/s	64–72 mph 56–63 knots 28.5–32.6 m/s	≥73 mph ≥63 knots ≥32.7 m/s

Vindens virkning på fjellet
Navn	m/s	knop	Kjennetegn
Stille	0,0-0,2	0-1	Snøfiller daler omtrent rett ned, gjerne i en pendlende bevegelse.
Flau vind	0,3-1,5	1-3	Så vidt følbar. Det er tydelig at snøfillene driver med vinden.
Svak vind	1,6-3,3	4-6	Godt følbar i sterk kulde. Snøfillene beveger seg mer horisontalt enn vertikalt.
Lett bris	3,4-5,4	7-10	Vinden merkes tydelig og kan sjenere. Fallende snø synes å bevege seg meget raskere horisontalt enn vertikalt.
Laber bris	5,5-7,9	11-16	Er ubehagelig i kaldt vær og gir merkbar motstand. Fallende snø hvirvler av sted med vinden. Snødrevet mot ansiktet er meget sjenerende.
Frisk bris	8,0-10,7	17-21	Det blir tungt å gå på ski mot været. Fokksnø som driver langs bakken, hvirvles så høyt at synsvidden nedsettes. Snødrevet pisker i ansiktet.
Liten kuling	10,8-13,8	22-27	Det er meget slitsomt å ta seg fram mot været. Snøfokk setter ned sikten til under 1 km. Vanskelig å holde ubeskyttet ansikt mot vinden i lengre tid. Folk flest bør ikke legge ut på tur over snaufjellet ved denne og høyere vindstyrker.
Stiv kuling	13,9-17,1	28-33	I motvind må en lute seg fram over skiene og legge stor kraft i stavtakene, selv på flat mark. Det kan være vanskelig å holde bena i vindrossene. Snøfokk setter ned sikten til få hundre meter. Vanskelig å orientere seg i terrenget. En skitur i fjellet ved vindstyrke 7 er en stor påkjenning for de fleste.
Sterk kuling	17,2-20,7	34-40	Fjellet står i kok. Kvister og lav fra trærne driver med vinden. Meget vanskelig å gå på skiene. Nesten umulig å bære skiene på nakken. Snøfokk setter ned sikten til under 100 m. Umulig å orientere seg i terrenget. Meget vanskelig å følge selv godt kvistede løyper. Legg ikke ut på tur!
Liten storm	20,8-24,4	41-47	Vind og snøfokk gjør det umulig å ta seg fram på ski over fjellet. Selv i klarvær og lite snøfokk kan påkjenningen bli så stor at en snøhule eller hytte er eneste redning.
Full storm	24,5-28,4	48-55	Denne og høyere vindstyrker vil de fleste aldri komme ut for. Trær velter over ledninger for telefon og strøm. Det knaker i tømmervegger. Lette småhus rives av grunnmuren.
Sterk storm	28,5-32,6	56-63	Veier og jernbanelinjer blokkeres. Det er kaos på telefon- og strømnettet. Skog blir rasert.
Orkan	32,6-	64-	Hvis bebyggelse rammes blir det en naturkatastrofe som gjerne vil kreve flere menneskeliv.

Imports



In [1]:

    
# -*- coding: utf-8 -*-
%matplotlib inline
# %matplotlib ipympl
import matplotlib.pyplot as plt
# plt.style.use(['dark_background'])
import matplotlib.patches as patches
# import pylab as plt
plt.rcParams['figure.figsize'] = (14, 6)
import seaborn as sns
# sns.set_style("darkgrid")

import datetime
import numpy as np
import netCDF4

from pathlib import Path
from windrose import WindroseAxes

import warnings
warnings.filterwarnings("ignore")









    



C:\ProgramData\Anaconda3\lib\site-packages\windrose\windrose.py:29: MatplotlibDeprecationWarning: 
The Appender class was deprecated in Matplotlib 3.1 and will be removed in 3.3.
  addendum = docstring.Appender(msg, "\n\n")
C:\ProgramData\Anaconda3\lib\site-packages\windrose\windrose.py:30: MatplotlibDeprecationWarning: 
The copy_dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use docstring.copy() and cbook.dedent() instead.
  return lambda func: addendum(docstring.copy_dedent(base)(func))
C:\ProgramData\Anaconda3\lib\site-packages\windrose\windrose.py:30: MatplotlibDeprecationWarning: 
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.getdoc() instead.
  return lambda func: addendum(docstring.copy_dedent(base)(func))
C:\ProgramData\Anaconda3\lib\site-packages\windrose\windrose.py:30: MatplotlibDeprecationWarning: 
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
  return lambda func: addendum(docstring.copy_dedent(base)(func))

Classification



In [2]:

    
# Classify by Beaufort scale
def classify_wind(wind_speed_f, wind_dir_f, js_str=True):
    len_wind = len(wind_speed_f)

    # Seperate by wind direction
    N = wind_speed_f[np.where((wind_dir_f>-22.5) & (wind_dir_f<=22.5))]

    NE = wind_speed_f[np.where((wind_dir_f>22.5) & (wind_dir_f<=67.5))]

    E = wind_speed_f[np.where((wind_dir_f>67.5) & (wind_dir_f<=112.5))]

    SE = wind_speed_f[np.where((wind_dir_f>112.5) & (wind_dir_f<=167.5))]

    S = wind_speed_f[np.where((wind_dir_f>167.5) & (wind_dir_f<=180.0) | (wind_dir_f>-167.5) & (wind_dir_f<=-180.0))]

    SW = wind_speed_f[np.where((wind_dir_f>-167.5) & (wind_dir_f<=-112.5))]

    W = wind_speed_f[np.where((wind_dir_f>-112.5) & (wind_dir_f<=-67.5))]

    NW = wind_speed_f[np.where((wind_dir_f>-67.5) & (wind_dir_f<=-22.5))]


    # Seperate by wind speed
    gentle_breeze = 5.5
    strong_breeze = 13.8
    gale = 20.7
    storm = 32.6

    N_gentle_breeze = len(N[np.where(N<gentle_breeze)]) / len_wind * 100
    N_strong_breeze = len(N[np.where((N>=gentle_breeze) & (N<strong_breeze))]) / len_wind * 100
    N_gale = len(N[np.where((N>=strong_breeze) & (N<gale))]) / len_wind * 100
    N_storm = len(N[np.where((N>=gale) & (N<storm))]) / len_wind * 100
    N_hurricane = len(N[np.where((N>=storm))]) / len_wind * 100

    NE_gentle_breeze = len(NE[np.where(NE<gentle_breeze)]) / len_wind * 100
    NE_strong_breeze = len(NE[np.where((NE>=gentle_breeze) & (NE<strong_breeze))]) / len_wind * 100
    NE_gale = len(NE[np.where((NE>=strong_breeze) & (NE<gale))]) / len_wind * 100
    NE_storm = len(NE[np.where((NE>=gale) & (NE<storm))]) / len_wind * 100
    NE_hurricane = len(NE[np.where((NE>=storm))]) / len_wind * 100

    E_gentle_breeze = len(E[np.where(E<gentle_breeze)]) / len_wind * 100
    E_strong_breeze = len(E[np.where((E>=gentle_breeze) & (E<strong_breeze))]) / len_wind * 100
    E_gale = len(E[np.where((E>=strong_breeze) & (E<gale))]) / len_wind * 100
    E_storm = len(E[np.where((E>=gale) & (E<storm))]) / len_wind * 100
    E_hurricane = len(E[np.where((E>=storm))]) / len_wind * 100

    SE_gentle_breeze = len(SE[np.where(SE<gentle_breeze)]) / len_wind * 100
    SE_strong_breeze = len(SE[np.where((SE>=gentle_breeze) & (SE<strong_breeze))]) / len_wind * 100
    SE_gale = len(SE[np.where((SE>=strong_breeze) & (SE<gale))]) / len_wind * 100
    SE_storm = len(SE[np.where((SE>=gale) & (SE<storm))]) / len_wind * 100
    SE_hurricane = len(SE[np.where((SE>=storm))]) / len_wind * 100

    S_gentle_breeze = len(S[np.where(S<gentle_breeze)]) / len_wind * 100
    S_strong_breeze = len(S[np.where((S>=gentle_breeze) & (S<strong_breeze))]) / len_wind * 100
    S_gale = len(S[np.where((S>=strong_breeze) & (S<gale))]) / len_wind * 100
    S_storm = len(S[np.where((S>=gale) & (S<storm))]) / len_wind * 100
    S_hurricane = len(S[np.where((S>=storm))]) / len_wind * 100

    SW_gentle_breeze = len(SW[np.where(SW<gentle_breeze)]) / len_wind * 100
    SW_strong_breeze = len(SW[np.where((SW>=gentle_breeze) & (SW<strong_breeze))]) / len_wind * 100
    SW_gale = len(SW[np.where((SW>=strong_breeze) & (SW<gale))]) / len_wind * 100
    SW_storm = len(SW[np.where((SW>=gale) & (SW<storm))]) / len_wind * 100
    SW_hurricane = len(SW[np.where((SW>=storm))]) / len_wind * 100

    W_gentle_breeze = len(W[np.where(W<gentle_breeze)]) / len_wind * 100
    W_strong_breeze = len(W[np.where((W>=gentle_breeze) & (W<strong_breeze))]) / len_wind * 100
    W_gale = len(W[np.where((W>=strong_breeze) & (W<gale))]) / len_wind * 100
    W_storm = len(W[np.where((W>=gale) & (W<storm))]) / len_wind * 100
    W_hurricane = len(W[np.where((W>=storm))]) / len_wind * 100

    NW_gentle_breeze = len(NW[np.where(NW<gentle_breeze)]) / len_wind * 100
    NW_strong_breeze = len(NW[np.where((NW>=gentle_breeze) & (NW<strong_breeze))]) / len_wind * 100
    NW_gale = len(NW[np.where((NW>=strong_breeze) & (NW<gale))]) / len_wind * 100
    NW_storm = len(NW[np.where((NW>=gale) & (NW<storm))]) / len_wind * 100
    NW_hurricane = len(NW[np.where((NW>=storm))]) / len_wind * 100

    #### can be removed after testing
    ####
    print([N_gentle_breeze, NE_gentle_breeze, E_gentle_breeze, SE_gentle_breeze, S_gentle_breeze, SW_gentle_breeze, W_gentle_breeze, NW_gentle_breeze])
    print([N_strong_breeze, NE_strong_breeze, E_strong_breeze, SE_strong_breeze, S_strong_breeze, SW_strong_breeze, W_strong_breeze, NW_strong_breeze])
    print([N_gale, NE_gale, E_gale, SE_gale, S_gale, SW_gale, W_gale, NW_gale])
    print([N_storm, NE_storm, E_storm, SE_storm, S_storm, SW_storm, W_storm, NW_storm])
    print([N_hurricane, NE_hurricane, E_hurricane, SE_hurricane, S_hurricane, SW_hurricane, W_hurricane, NW_hurricane])
    ####

    if js_str:
        str0 = "var light_winds = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_gentle_breeze, NE_gentle_breeze, E_gentle_breeze, SE_gentle_breeze, S_gentle_breeze, SW_gentle_breeze, W_gentle_breeze, NW_gentle_breeze)
        str1 = "var breeze = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_strong_breeze, NE_strong_breeze, E_strong_breeze, SE_strong_breeze, S_strong_breeze, SW_strong_breeze, W_strong_breeze, NW_strong_breeze)
        str2 = "var gale = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_gale, NE_gale, E_gale, SE_gale, S_gale, SW_gale, W_gale, NW_gale)
        str3 = "var storm = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_storm, NE_storm, E_storm, SE_storm, S_storm, SW_storm, W_storm, NW_storm)
        str4 = "var hurricane = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_hurricane, NE_hurricane, E_hurricane, SE_hurricane, S_hurricane, SW_hurricane, W_hurricane, NW_hurricane)

        json_str = "{0}\n{1}\n{2}\n{3}\n{4}".format(str0, str1, str2, str3, str4)

        return json_str

    else:
        wind_classes = {"gentle_breeze": [N_gentle_breeze, NE_gentle_breeze, E_gentle_breeze, SE_gentle_breeze, S_gentle_breeze, SW_gentle_breeze, W_gentle_breeze, NW_gentle_breeze],
                        "strong_breeze": [N_strong_breeze, NE_strong_breeze, E_strong_breeze, SE_strong_breeze, S_strong_breeze, SW_strong_breeze, W_strong_breeze, NW_strong_breeze],
                        "gale": [N_gale, NE_gale, E_gale, SE_gale, S_gale, SW_gale, W_gale, NW_gale],
                        "storm":[N_storm, NE_storm, E_storm, SE_storm, S_storm, SW_storm, W_storm, NW_storm],
                        "hurricane": [N_hurricane, NE_hurricane, E_hurricane, SE_hurricane, S_hurricane, SW_hurricane, W_hurricane, NW_hurricane]}

        return wind_classes



In [3]:

    
def classify_wind_transport(wind_speed_f, wind_dir_f, js_str=True):
    """
    Classes with regard to transport threshold.
    """
    len_wind = len(wind_speed_f)

    # Seperate by wind direction
    N = wind_speed_f[np.where((wind_dir_f>-22.5) & (wind_dir_f<=22.5))]

    NE = wind_speed_f[np.where((wind_dir_f>22.5) & (wind_dir_f<=67.5))]

    E = wind_speed_f[np.where((wind_dir_f>67.5) & (wind_dir_f<=112.5))]

    SE = wind_speed_f[np.where((wind_dir_f>112.5) & (wind_dir_f<=167.5))]

    S = wind_speed_f[np.where((wind_dir_f>167.5) & (wind_dir_f<=180.0) | (wind_dir_f>-167.5) & (wind_dir_f<=-180.0))]

    SW = wind_speed_f[np.where((wind_dir_f>-167.5) & (wind_dir_f<=-112.5))]

    W = wind_speed_f[np.where((wind_dir_f>-112.5) & (wind_dir_f<=-67.5))]

    NW = wind_speed_f[np.where((wind_dir_f>-67.5) & (wind_dir_f<=-22.5))]


    # Seperate by wind speed
    no_transport = 3.3 # less than light breeze, <2
    snowfall = 3.3 # more than light breeze, 3
    dry_snow = 5.5 # more than gentle breeze, 4
    wet_snow = 7.9 # more than moderate breeze, 5
    all_snow = 10.7 # more than fresh breeze, >5

    N_no_transport = len(N[np.where(N<no_transport)]) / len_wind * 100
    N_snowfall = len(N[np.where((N>=snowfall) & (N<dry_snow))]) / len_wind * 100
    N_dry_snow = len(N[np.where((N>=dry_snow) & (N<wet_snow))]) / len_wind * 100
    N_wet_snow = len(N[np.where((N>=wet_snow) & (N<all_snow))]) / len_wind * 100
    N_all_snow = len(N[np.where((N>=all_snow))]) / len_wind * 100

    NE_no_transport = len(NE[np.where(NE<no_transport)]) / len_wind * 100
    NE_snowfall = len(NE[np.where((NE>=snowfall) & (NE<dry_snow))]) / len_wind * 100
    NE_dry_snow = len(NE[np.where((NE>=dry_snow) & (NE<wet_snow))]) / len_wind * 100
    NE_wet_snow = len(NE[np.where((NE>=wet_snow) & (NE<all_snow))]) / len_wind * 100
    NE_all_snow = len(NE[np.where((NE>=all_snow))]) / len_wind * 100

    E_no_transport = len(E[np.where(E<no_transport)]) / len_wind * 100
    E_snowfall = len(E[np.where((E>=snowfall) & (E<dry_snow))]) / len_wind * 100
    E_dry_snow = len(E[np.where((E>=dry_snow) & (E<wet_snow))]) / len_wind * 100
    E_wet_snow = len(E[np.where((E>=wet_snow) & (E<all_snow))]) / len_wind * 100
    E_all_snow = len(E[np.where((E>=all_snow))]) / len_wind * 100

    SE_no_transport = len(SE[np.where(SE<no_transport)]) / len_wind * 100
    SE_snowfall = len(SE[np.where((SE>=snowfall) & (SE<dry_snow))]) / len_wind * 100
    SE_dry_snow = len(SE[np.where((SE>=dry_snow) & (SE<wet_snow))]) / len_wind * 100
    SE_wet_snow = len(SE[np.where((SE>=wet_snow) & (SE<all_snow))]) / len_wind * 100
    SE_all_snow = len(SE[np.where((SE>=all_snow))]) / len_wind * 100

    S_no_transport = len(S[np.where(S<no_transport)]) / len_wind * 100
    S_snowfall = len(S[np.where((S>=snowfall) & (S<dry_snow))]) / len_wind * 100
    S_dry_snow = len(S[np.where((S>=dry_snow) & (S<wet_snow))]) / len_wind * 100
    S_wet_snow = len(S[np.where((S>=wet_snow) & (S<all_snow))]) / len_wind * 100
    S_all_snow = len(S[np.where((S>=all_snow))]) / len_wind * 100

    SW_no_transport = len(SW[np.where(SW<no_transport)]) / len_wind * 100
    SW_snowfall = len(SW[np.where((SW>=snowfall) & (SW<dry_snow))]) / len_wind * 100
    SW_dry_snow = len(SW[np.where((SW>=dry_snow) & (SW<wet_snow))]) / len_wind * 100
    SW_wet_snow = len(SW[np.where((SW>=wet_snow) & (SW<all_snow))]) / len_wind * 100
    SW_all_snow = len(SW[np.where((SW>=all_snow))]) / len_wind * 100

    W_no_transport = len(W[np.where(W<no_transport)]) / len_wind * 100
    W_snowfall = len(W[np.where((W>=snowfall) & (W<dry_snow))]) / len_wind * 100
    W_dry_snow = len(W[np.where((W>=dry_snow) & (W<wet_snow))]) / len_wind * 100
    W_wet_snow = len(W[np.where((W>=wet_snow) & (W<all_snow))]) / len_wind * 100
    W_all_snow = len(W[np.where((W>=all_snow))]) / len_wind * 100

    NW_no_transport = len(NW[np.where(NW<no_transport)]) / len_wind * 100
    NW_snowfall = len(NW[np.where((NW>=snowfall) & (NW<dry_snow))]) / len_wind * 100
    NW_dry_snow = len(NW[np.where((NW>=dry_snow) & (NW<wet_snow))]) / len_wind * 100
    NW_wet_snow = len(NW[np.where((NW>=wet_snow) & (NW<all_snow))]) / len_wind * 100
    NW_all_snow = len(NW[np.where((NW>=all_snow))]) / len_wind * 100

    print([N_no_transport, NE_no_transport, E_no_transport, SE_no_transport, S_no_transport, SW_no_transport, W_no_transport, NW_no_transport])
    print([N_snowfall, NE_snowfall, E_snowfall, SE_snowfall, S_snowfall, SW_snowfall, W_snowfall, NW_snowfall])
    print([N_dry_snow, NE_dry_snow, E_dry_snow, SE_dry_snow, S_dry_snow, SW_dry_snow, W_dry_snow, NW_dry_snow])
    print([N_wet_snow, NE_wet_snow, E_wet_snow, SE_wet_snow, S_wet_snow, SW_wet_snow, W_wet_snow, NW_wet_snow])
    print([N_all_snow, NE_all_snow, E_all_snow, SE_all_snow, S_all_snow, SW_all_snow, W_all_snow, NW_all_snow])

    if js_str:
        str0 = "var no_transport = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_no_transport, NE_no_transport, E_no_transport, SE_no_transport, S_no_transport, SW_no_transport, W_no_transport, NW_no_transport)
        str1 = "var snowfall = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_snowfall, NE_snowfall, E_snowfall, SE_snowfall, S_snowfall, SW_snowfall, W_snowfall, NW_snowfall)
        str2 = "var dry_snow = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_dry_snow, NE_dry_snow, E_dry_snow, SE_dry_snow, S_dry_snow, SW_dry_snow, W_dry_snow, NW_dry_snow)
        str3 = "var wet_snow = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_wet_snow, NE_wet_snow, E_wet_snow, SE_wet_snow, S_wet_snow, SW_wet_snow, W_wet_snow, NW_wet_snow)
        str4 = "var all_snow = [{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}];".format(N_all_snow, NE_all_snow, E_all_snow, SE_all_snow, S_all_snow, SW_all_snow, W_all_snow, NW_all_snow)

        json_str = "{0}\n{1}\n{2}\n{3}\n{4}".format(str0, str1, str2, str3, str4)

        return json_str

    else:
        wind_classes = {"no_transport": [N_no_transport, NE_no_transport, E_no_transport, SE_no_transport, S_no_transport, SW_no_transport, W_no_transport, NW_no_transport],
                        "snowfall": [N_snowfall, NE_snowfall, E_snowfall, SE_snowfall, S_snowfall, SW_snowfall, W_snowfall, NW_snowfall],
                        "dry_snow": [N_dry_snow, NE_dry_snow, E_dry_snow, SE_dry_snow, S_dry_snow, SW_dry_snow, W_dry_snow, NW_dry_snow],
                        "wet_snow":[N_wet_snow, NE_wet_snow, E_wet_snow, SE_wet_snow, S_wet_snow, SW_wet_snow, W_wet_snow, NW_wet_snow],
                        "all_snow": [N_all_snow, NE_all_snow, E_all_snow, SE_all_snow, S_all_snow, SW_all_snow, W_all_snow, NW_all_snow]}

        return wind_classes



In [4]:

    
# Define plotting functions
def plt_wind_dir(wind_dir, ax=None):
    if ax == None:
        f, (ax) = plt.subplots(1, 1)
    pwd = ax.imshow(np.flipud(wind_dir), cmap=plt.cm.hsv, vmin=0, vmax=360)
    cbar = plt.colorbar(mappable=pwd, ax=ax)
    cbar.set_ticks([0, 45, 90, 135, 180, 225, 270, 315, 360])
    cbar.set_ticklabels(['N', 'NØ', 'Ø', 'SØ', 'S', 'SV', 'V', 'NV', 'N'])
    # cbar.set_ticks([-180, -135, -90, -45, 0, 45, 90, 135, 180])
    # cbar.set_ticklabels(['S', 'SV', 'V', 'NV', 'N', 'NØ', 'Ø', 'SØ', 'S'])
    cbar.set_label("Vindretning")

def plt_wind_speed(wind_speed, ax=None):
    if ax == None:
        f, (ax) = plt.subplots(1, 1)
    pws = ax.imshow(np.flipud(wind_speed), cmap=plt.cm.jet, vmin=0, vmax=32.6)
    cbar = plt.colorbar(mappable=pws, ax=ax)
    cbar.set_ticks([0, 1.5, 3.3, 5.5, 7.9, 10.7, 13.8, 17.1, 20.7, 24.4, 28.4, 32.6])
    cbar.set_label("Vindhastighet (m/s)")

def plt_windrose(wind_dir_f, wind_speed_f, title=""):
    # Convert wind direction from -180 to 180 to 0 to 360 degree.
    wind_dir_f[np.where((wind_dir_f<0.0))] = wind_dir_f[np.where((wind_dir_f<0.0))] + 360.0

    ax = WindroseAxes.from_ax()
    ax.bar(wind_dir_f, wind_speed_f, normed=True, opening=0.8, edgecolor='white', nsector=8, bins=[0,1.6,3.4,5.5,8.0,10.8,13.9,17.2,20.8,24.5,28.5,32.6])
    ax.set_legend()
    ax.set_title(title)
    ax.set_facecolor('k')
    plt.gcf().patch.set_facecolor('k')
    return ax

def plt_windrose_transport(wind_dir_f, wind_speed_f, title=""):
    # Convert wind direction from -180 to 180 to 0 to 360 degree.
    wind_dir_f[np.where((wind_dir_f<0.0))] = wind_dir_f[np.where((wind_dir_f<0.0))] + 360.0

    ax = WindroseAxes.from_ax()
    ax.bar(wind_dir_f, wind_speed_f, normed=True, opening=0.8, edgecolor='white', nsector=8, bins=[0,2,4.5,6.5,10.7,13.6])
    #ax.legend(labels=["ingen transport", "ved snøfall", "ved ubunden snø", "ved bunden snø", "hvis ikke skare/avblåst"])
    ax.set_legend()
    ax.set_title(title)
    ax.set_facecolor('k')
    return ax

def plt_snow_transport(wind_dir_f, wind_speed_f, title=""):
    # Convert wind direction from -180 to 180 to 0 to 360 degree.
    #wind_dir_f = -wind_dir_f
    wind_dir_f[np.where((wind_dir_f<0.0))] = wind_dir_f[np.where((wind_dir_f<0.0))] + 360.0

    ax = WindroseAxes.from_ax()
    ax.contourf(wind_dir_f, wind_speed_f, normed=True, bins=[0,2,4.5,6.5,10.7,13.6])
    ax.set_legend()
    ax.set_title(title)

    return ax

Short test



In [5]:

    
wind_dir = np.array([[0, 30, 50], [45, 180, 220], [350, 200, 130]], dtype="float")
print(wind_dir)
plt_wind_dir(wind_dir)









    



[[  0.  30.  50.]
 [ 45. 180. 220.]
 [350. 200. 130.]]

Using meps_det_pp files on



In [41]:

    
meps_dir = r'Y:\metdata\prognosis\met_pp_nordic\forecast\archive'
meps_yr = 2020
meps_file = 'met_forecast_1_0km_Norway_{date}T{hour}Z.nc'

meps_date = 20200416
meps_hour = 4

meps_pth = Path(meps_dir, str(meps_yr), meps_file.format(date=meps_date, hour=str(meps_hour).zfill(2)))
print(meps_pth)









    



Y:\metdata\prognosis\met_pp_nordic\forecast\archive\2020\met_forecast_1_0km_Norway_20200416T04Z.nc



In [42]:

    
nc = netCDF4.Dataset(meps_pth)

ref_time_v = nc.variables['forecast_reference_time']
ref_time = netCDF4.num2date(ref_time_v[0], ref_time_v.units)
print('Data from model run at {0}'.format(ref_time))

time_v = nc.variables['time']

print('Period from {0} to {1}'.format(netCDF4.num2date(time_v[0], time_v.units), netCDF4.num2date(time_v[-1], time_v.units)))

# Choose a time-step
t_index = 18
# Choose a pressure level (if applicable)
p_index = 12 # 12=1000hPa, 11=925hPa, 10=850hPa, ..., 7=500hPa, ..., 0=50hPa in arome_metcoop_test

ts = netCDF4.num2date(time_v[t_index], time_v.units)
print('Index chosen at time', ts)









    



Data from model run at 2020-04-16 04:00:00
Period from 2020-04-16 04:00:00 to 2020-04-18 13:00:00
Index chosen at time 2020-04-16 22:00:00



In [43]:

    
for k in nc.variables.keys():
    print(k)

wind_speed_v = nc.variables['wind_speed_10m']
wind_dir_v = nc.variables['wind_direction_10m']









    



forecast_reference_time
time
x
y
altitude
air_temperature_2m
precipitation_amount
wind_direction_10m
wind_speed_10m
cloud_area_fraction
air_pressure_at_sea_level
integral_of_surface_downwelling_shortwave_flux_in_air_wrt_time
relative_humidity_2m
projection_utm
lon
lat



In [ ]:

Wind speed



In [44]:

    
ws = wind_speed_v[t_index, :, :]
plt_wind_speed(ws)



In [45]:

    
sns.distplot(ws.flatten())









    Out[45]:





<matplotlib.axes._subplots.AxesSubplot at 0x25732681be0>

Wind direction



In [46]:

    
wd = wind_dir_v[t_index, :, :]
plt_wind_dir(wd)



In [47]:

    
sns.distplot(wd.flatten())









    Out[47]:





<matplotlib.axes._subplots.AxesSubplot at 0x25732767978>



In [48]:

    
# classify_wind(ws.flatten(), wd.flatten())



In [49]:

    
plt_windrose(wd.flatten(), ws.flatten(), str(ts))









    Out[49]:





<windrose.windrose.WindroseAxes at 0x257328ff9e8>



In [50]:

    
#plt_windrose_transport(wind_dir_f, wind_speed_f, t_info)



In [51]:

    
#plt_snow_transport(snow_dep_f, wind_speed_f)



In [52]:

    
# snow_classes = classify_wind_transport(wind_speed_f, snow_dep_f, js_str=True)
#print(snow_classes)

Using wind data on xgeo-grid

Possible issues:

I do not have to use the flipud function for the xgeo-gridded data. Does that influence the wind direction, i.e. is the y-vector in the opposite direction between the two data sets???



In [53]:

    
# Load region mask
vr = netCDF4.Dataset(r"../data/terrain_parameters/VarslingsOmr_2017.nc", "r")

regions = (vr.variables["VarslingsOmr_2017"][:]) #np.flipud - needs to be flipped up/down when using wind from netCDF files.

# ID = 3007 # Vest-Finnmark
ID = 3014 # Lofoten & Vesterålen
#ID = 3029 # Indre Sogn
#ID = 3033 # Hordalandskysten
region_mask = np.where(regions==ID)
# get the lower left and upper right corner of a rectangle around the region
y_min, y_max, x_min, x_max = min(region_mask[0].flatten()), max(region_mask[0].flatten()), min(region_mask[1].flatten()), max(region_mask[1].flatten())

plt.imshow(regions, vmin=3007, vmax=3040)









    Out[53]:





<matplotlib.image.AxesImage at 0x25732bd36d8>



In [54]:

    
x1, x2 = x_min, x_max # possible to add a buffer of step_x
y1, y2 = y_min, y_max # possible to add a buffer of step_y

print("Time period: {0} - {1}".format(netCDF4.num2date(time_v[t_index], time_v.units), netCDF4.num2date(time_v[t_index+24], time_v.units)))

# wind_x = wind_x_v[6:7, y1:y2, x1:x2]
# wind_y = wind_y_v[6:7, y1:y2, x1:x2]
# print(wind_x.shape)

wind_speed_reg = wind_speed_v[t_index:t_index+24, y1:y2, x1:x2]
wind_dir_reg = wind_dir_v[t_index:t_index+24, y1:y2, x1:x2]
print(wind_speed_reg.shape)

region_mask = regions[y1:y2, x1:x2] # redefine region_mask, now clipped to area of interest



wind_speed_reg_ma = np.ma.asarray(wind_speed_reg)
wind_dir_reg_ma = np.ma.asarray(wind_dir_reg)
for i in range(wind_speed_reg.shape[0]):
    wind_speed_reg_ma[i, :, :] = np.ma.masked_where(region_mask!=ID, wind_speed_reg_ma[i, :, :], copy=False)
    #print(type(wind_x[i, :, :]))
    wind_dir_reg_ma[i, :, :] = np.ma.masked_where(region_mask!=ID, wind_dir_reg_ma[i, :, :], copy=False)









    



Time period: 2020-04-16 22:00:00 - 2020-04-17 22:00:00
(24, 219, 199)



In [55]:

    
altitude_reg = nc.variables['altitude'][y1:y2, x1:x2]
sns.distplot(altitude_reg.flatten(), label='altitude')









    Out[55]:





<matplotlib.axes._subplots.AxesSubplot at 0x25732a249b0>



In [56]:

    
cutoff_elevation = 30.0
for i in range(wind_speed_reg.shape[0]):
    wind_speed_reg_ma[i, :, :] = np.ma.masked_where(altitude_reg<=cutoff_elevation, wind_speed_reg_ma[i, :, :], copy=False)
    #print(type(wind_x[i, :, :]))
    wind_dir_reg_ma[i, :, :] = np.ma.masked_where(altitude_reg<=cutoff_elevation, wind_dir_reg_ma[i, :, :], copy=False)

ToDo: split data in gridcells that are above and below treeline - use a fixed treeline elevation or the vegetation mask for each cell.



In [57]:

    
wind_dir_reg_f = wind_dir_reg_ma.flatten()
wind_speed_reg_f = wind_speed_reg_ma.flatten()



In [58]:

    
sns.distplot(wind_dir_reg_f)









    Out[58]:





<matplotlib.axes._subplots.AxesSubplot at 0x25732c6b8d0>



In [59]:

    
ax = sns.distplot(wind_speed_reg_f, label=str(ID))
ax.legend()









    Out[59]:





<matplotlib.legend.Legend at 0x257373df208>



In [60]:

    
# Choose a time-step
t_index_24h = 12
ts = netCDF4.num2date(time_v[t_index+t_index_24h], time_v.units)

f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
plt_wind_speed(wind_speed_reg_ma[t_index_24h,:,:], ax=ax1)
plt_wind_dir(wind_dir_reg_ma[t_index_24h,:,:], ax=ax2)

plt.title(ts)
plt.show()



In [61]:

    
plt_windrose(wind_dir_reg_f, wind_speed_reg_f)









    Out[61]:





<windrose.windrose.WindroseAxes at 0x25735cd54e0>

ToDo: Compare wind speeds for 10m and 850hpa - 850hpa might be too high to use as a basis!

NEXT

Define criteria for which wind speed and direction to choose for the regional weather.

Lower threshold for wind speeds of concern
Threshold in % when to consider only one wind direction
Which directions to include if several wind directions could be relevant?

CHECK Looks like there is a bug. Percentages do not add up to 100% and values are generally too low...

Wind direction for wind speeds above 7 m/s



In [62]:

    
dominant = wind_dir_reg_f[np.where(wind_speed_reg_f>7.0)]
sns.distplot(dominant)









    Out[62]:





<matplotlib.axes._subplots.AxesSubplot at 0x257348349e8>



In [28]:

    
# snow_classes = classify_wind_transport(wind_speed_f, snow_dep_f, js_str=True)
# print(snow_classes)

Testing influence of matrix convention on vector direction



In [29]:

    
def vec(x, y):
    vec_size = np.sqrt(x**2 + y**2)
    vec_dir = np.degrees(np.arctan2(x, y))

    return vec_size, vec_dir



In [30]:

    
x1 = np.array([[1, 2, -1],[-2, -1, 1],[2, -1, 1]])

y1 = np.array([[1, 2, -1],[-2, -1, 1],[2, -1, 1]])

f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
ax1.imshow(x1)
ax2.imshow(y1)









    Out[30]:





<matplotlib.image.AxesImage at 0x25735f11550>



In [31]:

    
f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
ax1.imshow(np.flipud(x1))
ax2.imshow(np.flipud(y1))









    Out[31]:





<matplotlib.image.AxesImage at 0x25731614128>



In [32]:

    
vec_size, vec_dir = vec(x1, y1)
plt_wind_dir(vec_dir)



In [33]:

    
vec_size, vec_dir = vec(np.flipud(x1), np.flipud(y1))
plt_wind_dir(vec_dir)



In [ ]: