WSL/SLF GitLab Repository

Commit 45a15df2 authored by Aschauer's avatar Aschauer
Browse files

add basic nan handling for `snowpack_evolution` function

parent fbb81d98
Pipeline #730 passed with stage
in 10 minutes and 6 seconds
...@@ -60,16 +60,16 @@ def _R_reduction_density( ...@@ -60,16 +60,16 @@ def _R_reduction_density(
): ):
R_r = np.zeros_like(rho_layers) R_r = np.zeros_like(rho_layers)
for i, r in enumerate(rho_layers): for i, r in enumerate(rho_layers):
if r <= rho_maxsettling: if r <= rho_maxsettling:
R_r[i] = 1 R_r[i] = 1
elif r >= rho_minsettling: elif r >= rho_minsettling:
R_r[i] = 0 R_r[i] = 0
else: else:
R_r[i] = (r - rho_minsettling)/(rho_maxsettling-rho_minsettling) R_r[i] = (r - rho_minsettling)/(rho_maxsettling-rho_minsettling)
return R_r return R_r
...@@ -87,13 +87,15 @@ def _get_settling_resistance( ...@@ -87,13 +87,15 @@ def _get_settling_resistance(
): ):
# calculate weights (scaling factors) of density and overburden. # calculate weights (scaling factors) of density and overburden.
w_rho = (R_max-R_min) * ratio_settling_influence_rho_vs_ob/(ratio_settling_influence_rho_vs_ob+1) w_rho = (R_max-R_min) * ratio_settling_influence_rho_vs_ob / \
(ratio_settling_influence_rho_vs_ob+1)
w_sigma = (R_max-R_min) * 1/(ratio_settling_influence_rho_vs_ob+1) w_sigma = (R_max-R_min) * 1/(ratio_settling_influence_rho_vs_ob+1)
# calculate reduction of R due to density and overburden. # calculate reduction of R due to density and overburden.
R_reduction_rho = _R_reduction_density(rho_layers, rho_minsettling, rho_maxsettling) R_reduction_rho = _R_reduction_density(rho_layers, rho_minsettling, rho_maxsettling)
R_reduction_sigma = _R_reduction_overburden(overburden_layers, ob_minsettling, ob_maxsettling) R_reduction_sigma = _R_reduction_overburden(
overburden_layers, ob_minsettling, ob_maxsettling)
R = R_max - w_rho*R_reduction_rho - w_sigma*R_reduction_sigma R = R_max - w_rho*R_reduction_rho - w_sigma*R_reduction_sigma
return R return R
...@@ -109,7 +111,7 @@ def _adjust_rho_max_wet_snowpack( ...@@ -109,7 +111,7 @@ def _adjust_rho_max_wet_snowpack(
swe_layers > 0, swe_layers > 0,
rho_max_wet - (rho_max_wet-rho_max_layers)*np.exp(-wetting_speed), rho_max_wet - (rho_max_wet-rho_max_layers)*np.exp(-wetting_speed),
rho_max_layers rho_max_layers
) )
@njit @njit
...@@ -146,20 +148,21 @@ def _adjust_rho_max_based_on_overburden( ...@@ -146,20 +148,21 @@ def _adjust_rho_max_based_on_overburden(
:class:`numpy.ndarray` :class:`numpy.ndarray`
Adapted rho_max in the layers based on overburden Adapted rho_max in the layers based on overburden
""" """
rho_max_current = np.zeros_like(overburden_layers) rho_max_current = np.zeros_like(overburden_layers)
for i, o in enumerate(overburden_layers): for i, o in enumerate(overburden_layers):
if o >= max_overburden: if o >= max_overburden:
rho_max_current[i] = rho_max_wet rho_max_current[i] = rho_max_wet
else: else:
rho_max_current[i] = ((rho_max_wet-rho_max_dry)/max_overburden)*o + rho_max_dry rho_max_current[i] = ((rho_max_wet-rho_max_dry) /
max_overburden)*o + rho_max_dry
updated_rho_max_layers = np.where((rho_max_current>rho_max_layers),
rho_max_current, updated_rho_max_layers = np.where((rho_max_current > rho_max_layers),
rho_max_current,
rho_max_layers) rho_max_layers)
updated_rho_max_layers = np.where((swe_layers>0), updated_rho_max_layers, 0.) updated_rho_max_layers = np.where((swe_layers > 0), updated_rho_max_layers, 0.)
return updated_rho_max_layers return updated_rho_max_layers
...@@ -185,11 +188,11 @@ def _remove_swe_from_top( ...@@ -185,11 +188,11 @@ def _remove_swe_from_top(
swe_removed = 0 swe_removed = 0
l = len(swe_layers)-1 l = len(swe_layers)-1
# melting from top means going backward in layers. # melting from top means going backward in layers.
while swe_removed > delta_swe: # both are (or will be) negative while swe_removed > delta_swe: # both are (or will be) negative
swe_removed = swe_removed - swe_layers[l] swe_removed = swe_removed - swe_layers[l]
swe_layers[l] = 0 swe_layers[l] = 0
l = l-1 l = l-1
# minimal floating point errors can cause the while loop to # minimal floating point errors can cause the while loop to
# run away, in that case we force it to stop at bottom of the # run away, in that case we force it to stop at bottom of the
# snowpack. # snowpack.
if l == -1: if l == -1:
...@@ -218,7 +221,7 @@ def timestep_forward( ...@@ -218,7 +221,7 @@ def timestep_forward(
overburden_minsettling, overburden_minsettling,
overburden_maxsettling, overburden_maxsettling,
max_overburden, max_overburden,
wetting_speed, wetting_speed,
): ):
"""Process a timestep forward. """Process a timestep forward.
...@@ -287,8 +290,8 @@ def timestep_forward( ...@@ -287,8 +290,8 @@ def timestep_forward(
rho_max_dry, rho_max_dry,
rho_max_wet, rho_max_wet,
max_overburden max_overburden
) )
# calculate settling resistance # calculate settling resistance
R = _get_settling_resistance( R = _get_settling_resistance(
rho_layers, rho_layers,
...@@ -300,14 +303,14 @@ def timestep_forward( ...@@ -300,14 +303,14 @@ def timestep_forward(
ob_maxsettling=overburden_maxsettling, ob_maxsettling=overburden_maxsettling,
rho_minsettling=rho_minsettling, rho_minsettling=rho_minsettling,
rho_maxsettling=rho_maxsettling, rho_maxsettling=rho_maxsettling,
) )
# update rho, i.e. calculate settling: # update rho, i.e. calculate settling:
rho_layers = np.where( rho_layers = np.where(
(swe_layers>0), (swe_layers > 0),
rho_max_layers - (rho_max_layers-rho_layers) * np.exp(-1/R), rho_max_layers - (rho_max_layers-rho_layers) * np.exp(-1/R),
0.) 0.)
if delta_swe > 0: if delta_swe > 0:
# rho_new should not get modified in the first timestep. # rho_new should not get modified in the first timestep.
rho_layers[-1] = rho_new rho_layers[-1] = rho_new
...@@ -315,6 +318,13 @@ def timestep_forward( ...@@ -315,6 +318,13 @@ def timestep_forward(
return swe_layers, rho_layers, rho_max_layers return swe_layers, rho_layers, rho_max_layers
@njit
def _set_layer_states_nan(swe_layers, rho_layers, rho_max_layers):
swe_layers = np.full_like(swe_layers, np.nan)
rho_layers = np.full_like(rho_layers, np.nan)
rho_max_layers = np.full_like(rho_max_layers, np.nan)
return swe_layers, rho_layers, rho_max_layers
@njit @njit
def _calculate_hs_layers(swe_layers, rho_layers): def _calculate_hs_layers(swe_layers, rho_layers):
...@@ -326,6 +336,29 @@ def _calculate_hs_layers(swe_layers, rho_layers): ...@@ -326,6 +336,29 @@ def _calculate_hs_layers(swe_layers, rho_layers):
return hs_layers return hs_layers
@njit
def _nansum_numba(array):
"""
Looped numba version faster than np.nansum
Parameters
----------
array : :class:`numpy.ndarray`
1 dimensional ndarray, float dtype
Returns
-------
sum : float
"""
sum = 0.
for x in array:
if np.isnan(x):
continue
else:
sum += x
return sum
@njit @njit
def _pad_end_of_array_with_zero(array): def _pad_end_of_array_with_zero(array):
# np.pad not supported in numba, we need some ugly hacking. # np.pad not supported in numba, we need some ugly hacking.
...@@ -365,7 +398,7 @@ def swe2hs_snowpack_evolution( ...@@ -365,7 +398,7 @@ def swe2hs_snowpack_evolution(
): ):
""" """
Snowpack evolution within the swe2hs model. Snowpack evolution within the swe2hs model.
Meant to be called on single hydrological years or chunks of nonzeros Meant to be called on single hydrological years or chunks of nonzeros
""" """
...@@ -380,25 +413,34 @@ def swe2hs_snowpack_evolution( ...@@ -380,25 +413,34 @@ def swe2hs_snowpack_evolution(
rho_max_layers_out = np.zeros((len(swe_input), len(swe_input))) rho_max_layers_out = np.zeros((len(swe_input), len(swe_input)))
# allocate layer containers # allocate layer containers
swe_layers = np.zeros(0) # tracking of swe swe_layers = np.zeros(0) # tracking of swe
rho_layers = np.zeros(0) # tracking of density rho_layers = np.zeros(0) # tracking of density
rho_max_layers = np.zeros(0) #tracking of maximum density rho_max_layers = np.zeros(0) # tracking of maximum density
# iterate through input array. # iterate through input array.
for i, swe in enumerate(swe_input): for i, swe in enumerate(swe_input):
# get delta swe:
if i==0:
delta_swe = swe
else:
delta_swe = swe - swe_input[i-1]
swe_layers, rho_layers, rho_max_layers = _pad_layer_arrays_with_zero( swe_layers, rho_layers, rho_max_layers = _pad_layer_arrays_with_zero(
swe_layers, swe_layers,
rho_layers, rho_layers,
rho_max_layers, rho_max_layers,
) )
if np.isnan(swe):
swe_layers, rho_layers, rho_max_layers = _set_layer_states_nan(
swe_layers,
rho_layers,
rho_max_layers
)
hs_out[i] = np.nan
continue
# get delta swe:
if i == 0:
delta_swe = swe
else:
delta_swe = swe - swe_input[i-1]
swe_layers, rho_layers, rho_max_layers = timestep_forward( swe_layers, rho_layers, rho_max_layers = timestep_forward(
delta_swe, delta_swe,
swe_layers, swe_layers,
...@@ -420,8 +462,8 @@ def swe2hs_snowpack_evolution( ...@@ -420,8 +462,8 @@ def swe2hs_snowpack_evolution(
hs_layers = _calculate_hs_layers(swe_layers, rho_layers) hs_layers = _calculate_hs_layers(swe_layers, rho_layers)
hs_out[i] = np.sum(hs_layers) hs_out[i] = _nansum_numba(hs_layers)
hs_layers_out[:i+1, i] = hs_layers hs_layers_out[:i+1, i] = hs_layers
rho_max_layers_out[:i+1, i] = rho_max_layers rho_max_layers_out[:i+1, i] = rho_max_layers
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment