| 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 |
In [1]:
# -*- coding: utf-8 -*-
%matplotlib inline
import matplotlib.patches as patches
import pylab as plt
plt.rcParams['figure.figsize'] = (14, 6)
import datetime
import numpy as np
import netCDF4
from windrose import WindroseAxes
import warnings
warnings.filterwarnings("ignore")
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=-180, vmax=180)
cbar = plt.colorbar(mappable=pwd, ax=ax)
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)
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)
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
In [5]:
wind_dir = np.degrees(np.arctan2(-np.array([[0, 20, 10], [-5, -18, -1], [1, -20, 3]], dtype="float"), -np.array([[10, 10, 0], [-5, 8, -10], [-11, 1, 3]], dtype="float")))
print(wind_dir)
plt_wind_dir(wind_dir)
In [9]:
#nc = netCDF4.Dataset("http://thredds.met.no/thredds/dodsC/arome25/arome_metcoop_default2_5km_latest.nc")
#nc = netCDF4.Dataset("http://thredds.met.no/thredds/dodsC/meps25files/meps_det_pp_2_5km_latest.nc")
#nc = netCDF4.Dataset("http://thredds.met.no/thredds/dodsC/meps25epsarchive/2017/08/13/meps_mbr0_pp_2_5km_20170813T00Z.nc")
nc = netCDF4.Dataset("http://thredds.met.no/thredds/dodsC/meps25epsarchive/2017/11/17/meps_mbr0_pp_2_5km_20171117T00Z.nc")
In [10]:
time_v = nc.variables['time']
# Choose a time-step
t_index = 6
# 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(ts)
In [11]:
wind_gust_v = nc.variables['wind_speed_of_gust']
plt_wind_speed(wind_gust_v[t_index, :, :])
In [12]:
wind_x_v = nc.variables['x_wind_10m']
wind_y_v = nc.variables['y_wind_10m']
In [13]:
f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
ax1.imshow(wind_x_v[t_index, :, :], cmap=plt.cm.seismic, vmin=-20, vmax=20)
ax2.imshow(wind_y_v[t_index, :, :], cmap=plt.cm.seismic, vmin=-20, vmax=20)
In [13]:
In [13]:
In [14]:
wind_speed = np.sqrt(wind_x_v[t_index, :, :]**2 + wind_y_v[t_index, :, :]**2)
plt_wind_speed(wind_speed)
The calculated vector indicates the direction to where the wind is blowing to, while we want to indicate the direction it is coming from in a wind map. Therefore we need to invert the x and y components of the vector before calculating the angle.
In [15]:
wind_dir = np.degrees(np.arctan2(-wind_x_v[t_index, :, :], -wind_y_v[t_index, :, :]))
plt_wind_dir(wind_dir)
In [15]:
In [16]:
### Lofoten
#x_range = np.arange(340, 430)
#y_range = np.arange(650,730)
### Langfjella
#x_range = np.arange(160, 260)
#y_range = np.arange(270,320)
### Hordalandskysten
x_range = np.arange(180, 220)
y_range = np.arange(290,370)
t_range = np.arange(18,42)
t1 = netCDF4.num2date(time_v[t_range[0]], time_v.units)
t2 = netCDF4.num2date(time_v[t_range[-1]], time_v.units)
t_info = "Period: {0} - {1}".format(t1, t2)
print(t_info)
plt.imshow((nc.variables["altitude"][:]))
rect = patches.Rectangle((x_range[0], y_range[0]), x_range[-1]-x_range[0], y_range[-1]-y_range[0], linewidth=1,edgecolor='r',facecolor='none')
plt.gca().add_patch(rect)
In [17]:
wind_speed = np.sqrt(wind_x_v[t_range, y_range, x_range]**2 + wind_y_v[t_range, y_range, x_range]**2)
wind_dir = np.degrees(np.arctan2(-wind_x_v[t_range, y_range, x_range], -wind_y_v[t_range, y_range, x_range]))
snow_dep = np.degrees(np.arctan2(wind_x_v[t_range, y_range, x_range], wind_y_v[t_range, y_range, x_range]))
In [18]:
laf = nc.variables["land_area_fraction"][y_range, x_range]
alt = nc.variables["altitude"][y_range, x_range]
f, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4)#, sharex=True, sharey=True)
ax1.imshow(np.flipud(laf))
ax2.imshow(np.flipud(alt))
plt_wind_speed(wind_speed[5,:,:], ax=ax3)
plt_wind_dir(wind_dir[5,:,:], ax=ax4)
In [19]:
print(wind_x_v[18, 300, 200], wind_y_v[18, 300, 200])
In [ ]:
In [19]:
In [20]:
f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharex=True, sharey=True)
ax1.imshow(np.flipud(wind_dir[5,:,:]), cmap=plt.cm.hsv,vmin=-180, vmax=180)
ax2.imshow(np.flipud(wind_dir[10,:,:]), cmap=plt.cm.hsv,vmin=-180, vmax=180)
ax3.imshow(np.flipud(wind_dir[15,:,:]), cmap=plt.cm.hsv,vmin=-180, vmax=180)
Out[20]:
In [21]:
wind_dir_f = wind_dir.flatten()
wind_speed_f = wind_speed.flatten()
snow_dep_f = snow_dep.flatten()
print(wind_dir_f.shape, wind_speed_f.shape)
len_wind = len(wind_speed_f)
In [22]:
classify_wind(wind_speed_f, wind_dir_f)
Out[22]:
In [23]:
plt_windrose(wind_dir_f, wind_speed_f, t_info)
Out[23]:
In [24]:
#plt_windrose_transport(wind_dir_f, wind_speed_f, t_info)
In [25]:
#plt_snow_transport(snow_dep_f, wind_speed_f)
In [26]:
snow_classes = classify_wind_transport(wind_speed_f, snow_dep_f, js_str=True)
#print(snow_classes)
In [23]:
# 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[23]:
In [24]:
#nc = netCDF4.Dataset(r"Y:\metdata\prognosis\arome\wind\netcdf\2017\arome2_5_wind_850_NVE_00_2017_04_27.nc", "r")
# nc = netCDF4.Dataset(r"Y:\metdata\prognosis\meps\det\archive\2017\mepsDet00_PTW_1km_20171117.nc", "r")
# nc = netCDF4.Dataset(r"Y:\metdata\prognosis\meps\det\archive\2020\meps_det_1km_20200319T03Z.nc", "r")
nc = netCDF4.Dataset(r"Y:\metdata\prognosis\meps\det\archive\2020\meps_det_1km_20200401T06Z.nc", "r")
time_v = nc.variables['time']
wind_x_v = nc.variables['x_wind_10m']
wind_y_v = nc.variables['y_wind_10m']
In [25]:
# Choose a time-step
t_index = 42
# Choose a pressure level (if applicable)
p_index = 12 # 12=1000hPa, 11=925hPa, 10=850hPa, ..., 7=500hPa, ..., 0=50hPa in arome_metcoop_test
print(netCDF4.num2date(time_v[t_index], time_v.units))
In [26]:
f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
ax1.imshow(wind_x_v[t_index, :, :], cmap=plt.cm.seismic, vmin=-20, vmax=20)
ax2.imshow(wind_y_v[t_index, :, :], cmap=plt.cm.seismic, vmin=-20, vmax=20)
Out[26]:
In [27]:
wind_speed = np.sqrt(wind_x_v[t_index, :, :]**2 + wind_y_v[t_index, :, :]**2)
plt_wind_speed(wind_speed)
In [28]:
wind_dir = np.degrees(np.arctan2(-wind_x_v[t_index, :, :], -wind_y_v[t_index, :, :]))
plt_wind_dir(wind_dir)
In [29]:
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_x = wind_x_v[t_index:t_index+24, y1:y2, x1:x2]
wind_y = wind_y_v[t_index:t_index+24, y1:y2, x1:x2]
print(wind_x.shape)
region_mask = regions[y1:y2, x1:x2] # redefine region_mask, now clipped to area of interest
wind_x_ma = np.ma.asarray(wind_x)
wind_y_ma = np.ma.asarray(wind_y)
for i in range(wind_x.shape[0]):
wind_x_ma[i, :, :] = np.ma.masked_where(region_mask!=ID, wind_x[i, :, :], copy=False)
#print(type(wind_x[i, :, :]))
wind_y_ma[i, :, :] = np.ma.masked_where(region_mask!=ID, wind_y[i, :, :], copy=False)
ToDo: split data in gridcells that are above and below treeline - use an fixed treeline elvation or the vegetation mask for each cell.
In [30]:
wind_speed_24h = np.sqrt(wind_x_ma**2 + wind_y_ma**2)
wind_dir_24h = np.degrees(np.arctan2(-wind_x_ma, -wind_y_ma)) # check if the minus is still applicable when using the xgeo-gridded data.
snow_dep_24h = np.degrees(np.arctan2(wind_x_ma, wind_y_ma)) # indicates the aspect where snow is deposited
print(type(wind_speed_24h))
wind_dir_f = wind_dir_24h.flatten()
wind_speed_f = wind_speed_24h.flatten()
snow_dep_f = snow_dep_24h.flatten()
print(wind_dir_f.shape, wind_speed_f.shape)
In [31]:
# Choose a time-step
t_index_24h = 1
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_24h[t_index_24h,:,:], ax=ax1)
plt_wind_dir(wind_dir_24h[t_index_24h,:,:], ax=ax2)
plt.title(ts)
plt.show()
In [32]:
plt_windrose(wind_dir_f, wind_speed_f)
Out[32]:
In [33]:
plt_windrose_transport(wind_dir_f, wind_speed_f)
Out[33]:
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.
CHECK Looks like there is a bug. Percentages do not add up to 100% and values are generally too low...
In [47]:
wind_classes = classify_wind_transport(wind_speed_f, wind_dir_f, js_str=True)
print(wind_classes)
In [48]:
dominant = wind_dir_f[np.where(wind_speed_f>7.0)]
plt.hist(dominant)
Out[48]:
In [49]:
snow_classes = classify_wind_transport(wind_speed_f, snow_dep_f, js_str=True)
print(snow_classes)
In [22]:
nc_url = "http://thredds.met.no/thredds/dodsC/arome25/arome_metcoop_test2_5km_latest.nc"
#nc_url = "http://thredds.met.no/thredds/dodsC/arome25/arome_metcoop_test2_5km_20170423_06.nc"
nc = netCDF4.Dataset(nc_url, 'r')
In [35]:
time_v = nc.variables['time']
pressure_v = nc.variables['pressure']
wind_x_v = nc.variables['x_wind_pl']
wind_y_v = nc.variables['y_wind_pl']
# Choose a time-period and spatial extent
x_range = np.arange(300, 400)
y_range = np.arange(600,750)
t_range = np.arange(6,30)
# Choose a pressure level (if applicable)
p_index = [12,11,10,7] # 12=1000hPa, 11=925hPa, 10=850hPa, ..., 7=500hPa, ..., 0=50hPa in arome_metcoop_test
In [36]:
t1 = netCDF4.num2date(time_v[t_range[0]], time_v.units)
t2 = netCDF4.num2date(time_v[t_range[-1]], time_v.units)
wind_speed = np.sqrt(wind_x_v[t_range, p_index, y_range, x_range]**2 + wind_y_v[t_range, p_index, y_range, x_range]**2)
wind_dir = np.degrees(np.arctan2(-wind_x_v[t_range, p_index, y_range, x_range], -wind_y_v[t_range, p_index, y_range, x_range]))
#f, axarr = plt.subplots(2, 2)
f = plt.figure(figsize=(12, 12))
i = 0
for pl in p_index: #p_index must be of length 4 at max
# Flatten arrays
wind_dir_f = wind_dir[:, i, :, :].flatten()
wind_speed_f = wind_speed[:, i, :, :].flatten()
print(wind_dir_f.shape, wind_speed_f.shape)
t_info = "Period: {0} - {1}\nPressure level: {2} hPa".format(t1, t2, pressure_v[pl])
ax = plt_windrose_transport(wind_dir_f, wind_speed_f, t_info)
#f.add_axes(ax)
#plt.savefig("windrose_pl_{0}.png".format(pl), dpi=90)
i+=1
TODO: Would be nice to get the plots equal in size.
In [42]:
plt_snow_transport(wind_dir_f, wind_speed_f)
Out[42]:
In [ ]:
In [ ]:
In [ ]:
In [136]:
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 [158]:
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[158]:
In [159]:
f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
ax1.imshow(np.flipud(x1))
ax2.imshow(np.flipud(y1))
Out[159]:
In [146]:
vec_size, vec_dir = vec(x1, y1)
plt_wind_dir(vec_dir)
In [144]:
vec_size, vec_dir = vec(np.flipud(x1), np.flipud(y1))
plt_wind_dir(vec_dir)
In [ ]: