From a7551194913a4e11fc2f80e5451e437bf7c6619f Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Wed, 29 Jul 2020 16:17:00 +0200 Subject: [PATCH 01/17] adding a messy code for iteration --- .../boundary_layer/RHmax_linearize_messy.py | 96 +++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py diff --git a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py new file mode 100644 index 0000000..a18736b --- /dev/null +++ b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py @@ -0,0 +1,96 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Jul 29 16:13:02 2020 + +Messy for now ... will have to be broken into individual parts as Leif as done with his refactoring +""" + +#%% +import os +import xarray as xr + +# load JOANNE dropsondes +input_dir = '/Users/annaleaalbright/Dropbox/EUREC4A/Dropsondes/Data/' +fp_dropsondes = os.path.join(input_dir,'EUREC4A_JOANNE_Dropsonde-RD41_Level_3_v0.5.7-alpha+0.g45fe69d.dirty.nc') +all_sondes = xr.open_dataset(fp_dropsondes) +#all_sondes = xr.open_dataset(fp_dropsondes).swap_dims({"sounding": "launch_time"}) + +#%% + +import numpy as np +import xarray as xr +from scipy.signal import find_peaks +import statsmodels.api as sm + +def calc_peakRH_linearize(ds, altitude, rh, time_dim, z_min, z_max, z_min_lin, z_max_lin): + """ + Calculate maximum in relative humidity (RH) profile that minimizes difference between + observed RH profile and RH profile linearized around lower levels (i.e. 200-400m). + Assume moisture is well-mixed at these lower levels, such that RH would increase and + assume linear profile as an idealization + Inputs: + -- ds: dataset + -- altitude, rh, time_dim: variable names + -- z_min and z_max: lower and upper bounds for mixed layer height + -- z_min_lin and z_max_lin: bounds for linearization of observed RH profile + Outputs: + -- da: datarray containing h_peakRH_linfit + """ + h_peakRH_linfit = np.zeros(len(ds[rh])) + + dz = int(ds[altitude].diff(dim=altitude)[1]) # dz=10m + + mixed_layer = np.logical_and(ds[altitude] >= z_min, ds[altitude] <= 1500) # enforce z_max later + ml_linfit = np.logical_and(ds[altitude] >= z_min_lin, ds[altitude] <= z_max_lin) # for linearized RH profile + + for i in range(len(ds[rh])): + + rh_profile = ds[rh].isel({time_dim: i}).interpolate_na(dim=altitude) + X_height = rh_profile[ml_linfit][altitude] + X = sm.add_constant(X_height.values) # add intercept + model = sm.OLS(rh_profile[ml_linfit].values,X).fit() # instantiate linear model + linearized_RH = model.predict(sm.add_constant(ds[altitude])) + + # 'find_peaks' is scipy function + idx_peaks_RH_raw, _ = find_peaks(rh_profile[mixed_layer]) + # add +X index from mixed layer (ml) lower bound (lb), divided by dz=10m + add_to_index_ml_lb = int(z_min / dz) + idx_peaks_RH = idx_peaks_RH_raw + add_to_index_ml_lb + obs_RH_peaks = rh_profile[idx_peaks_RH] # RH values at the local maxima + linearized_RH_values = linearized_RH[idx_peaks_RH] # linearized RH values at the observed RH maxima + min_diff_RH = min(abs(obs_RH_peaks - linearized_RH_values)) # find minimum RH difference between observed peaks and linearized profile + h_peakRH_linfit[i] = min_diff_RH[altitude] + + # set 'anomalously' high values > z_max to nan + # ... somewhat arbitrary choice + h_peakRH_linfit[h_peakRH_linfit > z_max] = np.nan + + dims = list(ds.dims.keys()) + # drop the height coord + del dims[dims.index(altitude)] + da = xr.DataArray(h_peakRH_linfit, dims=dims, coords={d: ds[d] for d in dims}) + da.attrs["long_name"] = "mixed layer height (from RH peak and linearization)" + da.attrs["units"] = "m" + return da + +#%% callfunction + +z_peakRH_linfit = calc_peakRH_linearize(all_sondes, altitude="height", rh="rh", time_dim="sounding", z_min=200.0, z_max=900.0, z_min_lin = 50.0, z_max_lin=200.0) + +#%% basic plots + +import seaborn as sns +import matplotlib.pyplot as plt +import numpy as np + +sns.set(context='notebook', style='whitegrid', palette='deep', font='sans-serif', font_scale=2, color_codes=True, rc=None) +plt.figure(figsize=(15,10)) +sns.distplot(z_peakRH_linfit) +plt.ylabel('pdf') +plt.xlabel('mixed layer height from RHmax linearize (m)') + +plt.figure(figsize=(15,10)) +plt.scatter(np.arange(len(z_peakRH_linfit)), z_peakRH_linfit) +plt.xlabel('sonde number') +plt.ylabel('mixed layer height from RHmax linearize (m)') \ No newline at end of file From edc6ebcad85f5a692463db5942968c4e46b655a6 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Wed, 29 Jul 2020 17:28:52 +0200 Subject: [PATCH 02/17] add messy file, BL height --- .../variables/boundary_layer/RHmax_linearize_messy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py index a18736b..392e3e0 100644 --- a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py +++ b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py @@ -23,7 +23,7 @@ from scipy.signal import find_peaks import statsmodels.api as sm -def calc_peakRH_linearize(ds, altitude, rh, time_dim, z_min, z_max, z_min_lin, z_max_lin): +def calc_peakRH_linearize(ds, altitude="height", rh="rh", time_dim="sounding", z_min, z_max, z_min_lin, z_max_lin): """ Calculate maximum in relative humidity (RH) profile that minimizes difference between observed RH profile and RH profile linearized around lower levels (i.e. 200-400m). From 9bb5db25d9cbd910c494cd399e035ea9479af47a Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Wed, 29 Jul 2020 17:39:23 +0200 Subject: [PATCH 03/17] minor change --- .gitignore | 1 + .../variables/boundary_layer/RHmax_linearize_messy.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 6e0bb90..509252a 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ .DS_Store *.nc data/ +eurec4a_environment.egg-info/ diff --git a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py index 392e3e0..79a560f 100644 --- a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py +++ b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py @@ -23,7 +23,7 @@ from scipy.signal import find_peaks import statsmodels.api as sm -def calc_peakRH_linearize(ds, altitude="height", rh="rh", time_dim="sounding", z_min, z_max, z_min_lin, z_max_lin): +def calc_peakRH_linearize(ds, altitude="height", rh="rh", time_dim="sounding", z_min=200.0, z_max=900.0, z_min_lin=50.0, z_max_lin=400.0): """ Calculate maximum in relative humidity (RH) profile that minimizes difference between observed RH profile and RH profile linearized around lower levels (i.e. 200-400m). From eec0c98c681f6d43f8693637769ed1edc95ff74c Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Wed, 29 Jul 2020 17:43:43 +0200 Subject: [PATCH 04/17] minor edit --- .../boundary_layer/RHmax_linearize_messy.py | 142 +++++++++++------- 1 file changed, 91 insertions(+), 51 deletions(-) diff --git a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py index 79a560f..6b7a2f0 100644 --- a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py +++ b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py @@ -11,10 +11,12 @@ import xarray as xr # load JOANNE dropsondes -input_dir = '/Users/annaleaalbright/Dropbox/EUREC4A/Dropsondes/Data/' -fp_dropsondes = os.path.join(input_dir,'EUREC4A_JOANNE_Dropsonde-RD41_Level_3_v0.5.7-alpha+0.g45fe69d.dirty.nc') +input_dir = "/Users/annaleaalbright/Dropbox/EUREC4A/Dropsondes/Data/" +fp_dropsondes = os.path.join( + input_dir, "EUREC4A_JOANNE_Dropsonde-RD41_Level_3_v0.5.7-alpha+0.g45fe69d.dirty.nc" +) all_sondes = xr.open_dataset(fp_dropsondes) -#all_sondes = xr.open_dataset(fp_dropsondes).swap_dims({"sounding": "launch_time"}) +# all_sondes = xr.open_dataset(fp_dropsondes).swap_dims({"sounding": "launch_time"}) #%% @@ -23,8 +25,18 @@ from scipy.signal import find_peaks import statsmodels.api as sm -def calc_peakRH_linearize(ds, altitude="height", rh="rh", time_dim="sounding", z_min=200.0, z_max=900.0, z_min_lin=50.0, z_max_lin=400.0): - """ + +def calc_peakRH_linearize( + ds, + altitude="height", + rh="rh", + time_dim="sounding", + z_min=200.0, + z_max=900.0, + z_min_lin=50.0, + z_max_lin=400.0, +): + """ Calculate maximum in relative humidity (RH) profile that minimizes difference between observed RH profile and RH profile linearized around lower levels (i.e. 200-400m). Assume moisture is well-mixed at these lower levels, such that RH would increase and @@ -37,60 +49,88 @@ def calc_peakRH_linearize(ds, altitude="height", rh="rh", time_dim="sounding", z Outputs: -- da: datarray containing h_peakRH_linfit """ - h_peakRH_linfit = np.zeros(len(ds[rh])) - - dz = int(ds[altitude].diff(dim=altitude)[1]) # dz=10m - - mixed_layer = np.logical_and(ds[altitude] >= z_min, ds[altitude] <= 1500) # enforce z_max later - ml_linfit = np.logical_and(ds[altitude] >= z_min_lin, ds[altitude] <= z_max_lin) # for linearized RH profile - - for i in range(len(ds[rh])): - - rh_profile = ds[rh].isel({time_dim: i}).interpolate_na(dim=altitude) - X_height = rh_profile[ml_linfit][altitude] - X = sm.add_constant(X_height.values) # add intercept - model = sm.OLS(rh_profile[ml_linfit].values,X).fit() # instantiate linear model - linearized_RH = model.predict(sm.add_constant(ds[altitude])) - - # 'find_peaks' is scipy function - idx_peaks_RH_raw, _ = find_peaks(rh_profile[mixed_layer]) - # add +X index from mixed layer (ml) lower bound (lb), divided by dz=10m - add_to_index_ml_lb = int(z_min / dz) - idx_peaks_RH = idx_peaks_RH_raw + add_to_index_ml_lb - obs_RH_peaks = rh_profile[idx_peaks_RH] # RH values at the local maxima - linearized_RH_values = linearized_RH[idx_peaks_RH] # linearized RH values at the observed RH maxima - min_diff_RH = min(abs(obs_RH_peaks - linearized_RH_values)) # find minimum RH difference between observed peaks and linearized profile - h_peakRH_linfit[i] = min_diff_RH[altitude] - - # set 'anomalously' high values > z_max to nan - # ... somewhat arbitrary choice - h_peakRH_linfit[h_peakRH_linfit > z_max] = np.nan - - dims = list(ds.dims.keys()) - # drop the height coord - del dims[dims.index(altitude)] - da = xr.DataArray(h_peakRH_linfit, dims=dims, coords={d: ds[d] for d in dims}) - da.attrs["long_name"] = "mixed layer height (from RH peak and linearization)" - da.attrs["units"] = "m" - return da + h_peakRH_linfit = np.zeros(len(ds[rh])) + + dz = int(ds[altitude].diff(dim=altitude)[1]) # dz=10m + + mixed_layer = np.logical_and( + ds[altitude] >= z_min, ds[altitude] <= 1500 + ) # enforce z_max later + ml_linfit = np.logical_and( + ds[altitude] >= z_min_lin, ds[altitude] <= z_max_lin + ) # for linearized RH profile + + for i in range(len(ds[rh])): + + rh_profile = ds[rh].isel({time_dim: i}).interpolate_na(dim=altitude) + X_height = rh_profile[ml_linfit][altitude] + X = sm.add_constant(X_height.values) # add intercept + model = sm.OLS( + rh_profile[ml_linfit].values, X + ).fit() # instantiate linear model + linearized_RH = model.predict(sm.add_constant(ds[altitude])) + + # 'find_peaks' is scipy function + idx_peaks_RH_raw, _ = find_peaks(rh_profile[mixed_layer]) + # add +X index from mixed layer (ml) lower bound (lb), divided by dz=10m + add_to_index_ml_lb = int(z_min / dz) + idx_peaks_RH = idx_peaks_RH_raw + add_to_index_ml_lb + obs_RH_peaks = rh_profile[idx_peaks_RH] # RH values at the local maxima + linearized_RH_values = linearized_RH[ + idx_peaks_RH + ] # linearized RH values at the observed RH maxima + min_diff_RH = min( + abs(obs_RH_peaks - linearized_RH_values) + ) # find minimum RH difference between observed peaks and linearized profile + h_peakRH_linfit[i] = min_diff_RH[altitude] + + # set 'anomalously' high values > z_max to nan + # ... somewhat arbitrary choice + h_peakRH_linfit[h_peakRH_linfit > z_max] = np.nan + + dims = list(ds.dims.keys()) + # drop the height coord + del dims[dims.index(altitude)] + da = xr.DataArray(h_peakRH_linfit, dims=dims, coords={d: ds[d] for d in dims}) + da.attrs["long_name"] = "mixed layer height (from RH peak and linearization)" + da.attrs["units"] = "m" + return da + #%% callfunction - -z_peakRH_linfit = calc_peakRH_linearize(all_sondes, altitude="height", rh="rh", time_dim="sounding", z_min=200.0, z_max=900.0, z_min_lin = 50.0, z_max_lin=200.0) + +z_peakRH_linfit = calc_peakRH_linearize( + all_sondes, + altitude="height", + rh="rh", + time_dim="sounding", + z_min=200.0, + z_max=900.0, + z_min_lin=50.0, + z_max_lin=200.0, +) #%% basic plots -import seaborn as sns +import seaborn as sns import matplotlib.pyplot as plt import numpy as np -sns.set(context='notebook', style='whitegrid', palette='deep', font='sans-serif', font_scale=2, color_codes=True, rc=None) -plt.figure(figsize=(15,10)) +sns.set( + context="notebook", + style="whitegrid", + palette="deep", + font="sans-serif", + font_scale=2, + color_codes=True, + rc=None, +) +plt.figure(figsize=(15, 10)) sns.distplot(z_peakRH_linfit) -plt.ylabel('pdf') -plt.xlabel('mixed layer height from RHmax linearize (m)') +plt.ylabel("pdf") +plt.xlabel("mixed layer height from RHmax linearize (m)") -plt.figure(figsize=(15,10)) +plt.figure(figsize=(15, 10)) plt.scatter(np.arange(len(z_peakRH_linfit)), z_peakRH_linfit) -plt.xlabel('sonde number') -plt.ylabel('mixed layer height from RHmax linearize (m)') \ No newline at end of file +plt.xlabel("sonde number") +plt.ylabel("mixed layer height from RHmax linearize (m)") From 3322c8666f42d7c9a072ee8ecd3c7c0b2aca343c Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Thu, 30 Jul 2020 09:35:09 +0200 Subject: [PATCH 05/17] Update calculate_density.py --- eurec4a_environment/variables/calculate_density.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/eurec4a_environment/variables/calculate_density.py b/eurec4a_environment/variables/calculate_density.py index 41081f4..462c091 100644 --- a/eurec4a_environment/variables/calculate_density.py +++ b/eurec4a_environment/variables/calculate_density.py @@ -1,18 +1,5 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -""" -Created on Thu Jul 30 09:29:04 2020 - -@author: annaleaalbright -""" - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Thu Jul 30 08:57:37 2020 - -@author: annaleaalbright -""" import numpy as np import xarray as xr From 686adf11627607e24de170c09cd9ec1ddadcd553 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Thu, 30 Jul 2020 15:36:13 +0200 Subject: [PATCH 06/17] test for calc_density function --- tests/test_density.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 tests/test_density.py diff --git a/tests/test_density.py b/tests/test_density.py new file mode 100644 index 0000000..4bd9515 --- /dev/null +++ b/tests/test_density.py @@ -0,0 +1,16 @@ +# -*- coding: utf-8 -*- + +import numpy as np + +from eurec4a_environment.variables.calc_density import calc_density +import eurec4a_environment.source_data + + +def test_density_calculation(): + ds = eurec4a_environment.source_data.open_joanne_dataset() + ds = ds.isel(sounding=slice(0, 10)) + da_density = calc_density(ds) + # check that we get some sensible numbers out for the density + assert da_density.units == "kg/m3" + assert da_density.max() < 2.0 + assert da_density.min() > 0.0 From f2312e6ca24ebfc4528468e9a75b74d798fd0b60 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Thu, 30 Jul 2020 15:38:41 +0200 Subject: [PATCH 07/17] fixed coordinates --- eurec4a_environment/variables/calculate_density.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/eurec4a_environment/variables/calculate_density.py b/eurec4a_environment/variables/calculate_density.py index 462c091..ba24075 100644 --- a/eurec4a_environment/variables/calculate_density.py +++ b/eurec4a_environment/variables/calculate_density.py @@ -3,10 +3,10 @@ import numpy as np import xarray as xr -from ...constants import Rd, Rv, eps +#from ...constants import Rd, Rv, eps +from constants import Rd, Rv, eps - -def calc_density(ds, pres="p", temp="T", specific_humidity="q"): +def calc_density(ds, pres="p", temp="T", specific_humidity="q", altitude="height"): # equation: rho = P/(Rd * Tv), where Tv = T(1 + mr/eps)/(1+mr) @@ -31,9 +31,11 @@ def calc_density(ds, pres="p", temp="T", specific_humidity="q"): density = (pressure) / ( Rd * (temp_K) * (1 + (mixing_ratio / eps)) / (1 + mixing_ratio) ) + density = density.transpose(transpose_coords=True) dims = list(ds.dims.keys()) da = xr.DataArray(density, dims=dims, coords={d: ds[d] for d in dims}) da.attrs["long_name"] = "density of air" da.attrs["units"] = "kg/m3" + return da From 22ba43ed98d246731209892b76c80230af326bf0 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Thu, 30 Jul 2020 15:53:03 +0200 Subject: [PATCH 08/17] minor edit --- .../variables/boundary_layer/RHmax_linearize_messy.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py index 6b7a2f0..bb616c2 100644 --- a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py +++ b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py @@ -3,20 +3,19 @@ """ Created on Wed Jul 29 16:13:02 2020 -Messy for now ... will have to be broken into individual parts as Leif as done with his refactoring """ #%% -import os -import xarray as xr - +#import os +#import xarray as xr +# # load JOANNE dropsondes input_dir = "/Users/annaleaalbright/Dropbox/EUREC4A/Dropsondes/Data/" fp_dropsondes = os.path.join( input_dir, "EUREC4A_JOANNE_Dropsonde-RD41_Level_3_v0.5.7-alpha+0.g45fe69d.dirty.nc" ) all_sondes = xr.open_dataset(fp_dropsondes) -# all_sondes = xr.open_dataset(fp_dropsondes).swap_dims({"sounding": "launch_time"}) +all_sondes = xr.open_dataset(fp_dropsondes).swap_dims({"sounding": "launch_time"}) #%% From 364a2c4744c81bfeefcc56192df582df67402be0 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Thu, 30 Jul 2020 16:04:48 +0200 Subject: [PATCH 09/17] add test file --- eurec4a_environment/variables/calculate_density.py | 1 + tests/test_density.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/eurec4a_environment/variables/calculate_density.py b/eurec4a_environment/variables/calculate_density.py index ba24075..4e146f2 100644 --- a/eurec4a_environment/variables/calculate_density.py +++ b/eurec4a_environment/variables/calculate_density.py @@ -4,6 +4,7 @@ import numpy as np import xarray as xr #from ...constants import Rd, Rv, eps + from constants import Rd, Rv, eps def calc_density(ds, pres="p", temp="T", specific_humidity="q", altitude="height"): diff --git a/tests/test_density.py b/tests/test_density.py index 4bd9515..c23237b 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -12,5 +12,5 @@ def test_density_calculation(): da_density = calc_density(ds) # check that we get some sensible numbers out for the density assert da_density.units == "kg/m3" - assert da_density.max() < 2.0 + assert da_density.max() < 2.5 assert da_density.min() > 0.0 From f7e1075855f115d8a3c18f0472216468a6405f15 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Thu, 30 Jul 2020 16:58:12 +0200 Subject: [PATCH 10/17] updating files --- .../boundary_layer/RHmax_linearize_messy.py | 4 +- .../boundary_layer/mixed_layer_height.py | 83 ++++++++++++++++- .../mixed_layer_height_from_gradient.py | 88 +++++++++++++++++++ .../{calculate_density.py => calc_density.py} | 5 +- 4 files changed, 172 insertions(+), 8 deletions(-) create mode 100644 eurec4a_environment/variables/boundary_layer/mixed_layer_height_from_gradient.py rename eurec4a_environment/variables/{calculate_density.py => calc_density.py} (96%) diff --git a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py index bb616c2..87d7702 100644 --- a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py +++ b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py @@ -6,8 +6,8 @@ """ #%% -#import os -#import xarray as xr +# import os +# import xarray as xr # # load JOANNE dropsondes input_dir = "/Users/annaleaalbright/Dropbox/EUREC4A/Dropsondes/Data/" diff --git a/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py b/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py index 60246ae..6c8a13c 100644 --- a/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py +++ b/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py @@ -30,7 +30,7 @@ def calc_peakRH_linearize( rh="rh", time_dim="sounding", z_min=200.0, - z_max=900.0, + z_max=1500.0, z_min_lin=50.0, z_max_lin=400.0, ): @@ -51,9 +51,7 @@ def calc_peakRH_linearize( dz = int(ds[altitude].diff(dim=altitude)[1]) # dz=10m - mixed_layer = np.logical_and( - ds[altitude] >= z_min, ds[altitude] <= 1500 - ) # enforce z_max later + mixed_layer = np.logical_and(ds[altitude] >= z_min, ds[altitude] <= z_max) ml_linfit = np.logical_and( ds[altitude] >= z_min_lin, ds[altitude] <= z_max_lin ) # for linearized RH profile @@ -93,3 +91,80 @@ def calc_peakRH_linearize( da.attrs["long_name"] = "mixed layer height (from RH peak and linearization)" da.attrs["units"] = "m" return da + + +def calc_from_gradient( + ds, var, threshold, z_min=200, altitude="height", time="sounding" +): + """ + Find mixed layer height as layer over which x(z+1) - x(z) < threshold + + Inputs: + -- ds: Dataset + -- var: variable name + -- threshold: in units of variable (i.e. g/kg or K). Maximum difference allowed in one vertical step, dz=10m, i.e. threshold_q = 0.4 g/kg, + or vertical gradient dq/dz < 0.04 g/kg/m + + Outputs: DataArray containing mixed layer height from gradient method + + Note that function is quite slow and should be optimized + + """ + + def calculateHmix_var( + density_profile, + var_profile, + threshold, + z_min, + altitude="height", + time="sounding", + ): + + var_diff = 0 + numer = 0 + denom = 0 + var_mix = 0 + + k = int(z_min / 10) # enforce lower limit = 200m. k is index of lower limit + + # var(z) - weighted_mean(z-1) < threshold + while abs(var_diff) < threshold: + + numer += 0.5 * ( + density_profile[k + 1] * var_profile[k + 1] + + density_profile[k] * var_profile[k] + ) + denom += 0.5 * (density_profile[k + 1] + density_profile[k]) + var_mix = numer / denom + k += 1 + var_diff = (var_profile[k] - var_mix).values + + hmix = var_profile[altitude].values[k] + return hmix + + # call inner function + from eurec4a_environment.variables.calc_density import calc_density + + da_density = calc_density(ds) + + hmix_vec = np.zeros(len(ds[var])) + for i in range((len(ds[var]))): + density_profile = da_density.isel({time: i}) + var_profile = ds[var].isel({time: i}) + hmix_vec[i] = calculateHmix_var( + density_profile, + var_profile, + threshold, + z_min, + altitude="height", + time="sounding", + ) + + dims = list(ds.dims.keys()) + # drop the height coord + del dims[dims.index(altitude)] + da = xr.DataArray(hmix_vec, dims=dims, coords={d: ds[d] for d in dims}) + da.attrs["long_name"] = f"mixed layer height (from gradient method using {var})" + da.attrs["units"] = "m" + + return da diff --git a/eurec4a_environment/variables/boundary_layer/mixed_layer_height_from_gradient.py b/eurec4a_environment/variables/boundary_layer/mixed_layer_height_from_gradient.py new file mode 100644 index 0000000..49868a9 --- /dev/null +++ b/eurec4a_environment/variables/boundary_layer/mixed_layer_height_from_gradient.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Find mixed layer height as layer over which x(z+1) - x(z) < threshold + +Inputs: + -- ds: Dataset + -- var: variable name + -- threshold: in units of variable (i.e. g/kg or K). Maximum difference allowed in one vertical step, dz=10m, i.e. threshold_q = 0.4 g/kg, + or vertical gradient dq/dz < 0.04 g/kg/m + +Outputs: DataArray containing mixed layer height from gradient method + +Note that function is quite slow and should be optimized + +""" + +import numpy as np +import xarray as xr +from scipy.signal import find_peaks +import statsmodels.api as sm +from constants import Rd, Rv, eps +import seaborn as sns + +# import eurec4a_environment.source_data +# ds = eurec4a_environment.source_data.open_joanne_dataset() + + +def calc_from_gradient( + ds, var, threshold, z_min=200, altitude="height", time="sounding" +): + def calculateHmix_var( + density_profile, + var_profile, + threshold, + z_min, + altitude="height", + time="sounding", + ): + + var_diff = 0 + numer = 0 + denom = 0 + var_mix = 0 + + k = int(z_min / 10) # enforce lower limit = 200m. k is index of lower limit + + # var(z) - weighted_mean(z-1) < threshold + while abs(var_diff) < threshold: + + numer += 0.5 * ( + density_profile[k + 1] * var_profile[k + 1] + + density_profile[k] * var_profile[k] + ) + denom += 0.5 * (density_profile[k + 1] + density_profile[k]) + var_mix = numer / denom + k += 1 + var_diff = (var_profile[k] - var_mix).values + + hmix = var_profile[altitude].values[k] + return hmix + + # call inner function + from eurec4a_environment.variables.calc_density import calc_density + + da_density = calc_density(ds) + + hmix_vec = np.zeros(len(ds[var])) + for i in range((len(ds[var]))): + density_profile = da_density.isel({time: i}) + var_profile = ds[var].isel({time: i}) + hmix_vec[i] = calculateHmix_var( + density_profile, + var_profile, + threshold, + z_min, + altitude="height", + time="sounding", + ) + + dims = list(ds.dims.keys()) + # drop the height coord + del dims[dims.index(altitude)] + da = xr.DataArray(hmix_vec, dims=dims, coords={d: ds[d] for d in dims}) + da.attrs["long_name"] = f"mixed layer height (from gradient method using {var})" + da.attrs["units"] = "m" + + return da diff --git a/eurec4a_environment/variables/calculate_density.py b/eurec4a_environment/variables/calc_density.py similarity index 96% rename from eurec4a_environment/variables/calculate_density.py rename to eurec4a_environment/variables/calc_density.py index 4e146f2..bbda8ca 100644 --- a/eurec4a_environment/variables/calculate_density.py +++ b/eurec4a_environment/variables/calc_density.py @@ -3,10 +3,11 @@ import numpy as np import xarray as xr -#from ...constants import Rd, Rv, eps +# from ...constants import Rd, Rv, eps from constants import Rd, Rv, eps + def calc_density(ds, pres="p", temp="T", specific_humidity="q", altitude="height"): # equation: rho = P/(Rd * Tv), where Tv = T(1 + mr/eps)/(1+mr) @@ -38,5 +39,5 @@ def calc_density(ds, pres="p", temp="T", specific_humidity="q", altitude="height da = xr.DataArray(density, dims=dims, coords={d: ds[d] for d in dims}) da.attrs["long_name"] = "density of air" da.attrs["units"] = "kg/m3" - + return da From 675a4b26905509983c469dd03789d47e1590b12f Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Thu, 30 Jul 2020 17:13:10 +0200 Subject: [PATCH 11/17] remove redundant code --- .../boundary_layer/RHmax_linearize_messy.py | 135 ------------------ eurec4a_environment/variables/calc_density.py | 5 +- 2 files changed, 2 insertions(+), 138 deletions(-) delete mode 100644 eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py diff --git a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py b/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py deleted file mode 100644 index 87d7702..0000000 --- a/eurec4a_environment/variables/boundary_layer/RHmax_linearize_messy.py +++ /dev/null @@ -1,135 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed Jul 29 16:13:02 2020 - -""" - -#%% -# import os -# import xarray as xr -# -# load JOANNE dropsondes -input_dir = "/Users/annaleaalbright/Dropbox/EUREC4A/Dropsondes/Data/" -fp_dropsondes = os.path.join( - input_dir, "EUREC4A_JOANNE_Dropsonde-RD41_Level_3_v0.5.7-alpha+0.g45fe69d.dirty.nc" -) -all_sondes = xr.open_dataset(fp_dropsondes) -all_sondes = xr.open_dataset(fp_dropsondes).swap_dims({"sounding": "launch_time"}) - -#%% - -import numpy as np -import xarray as xr -from scipy.signal import find_peaks -import statsmodels.api as sm - - -def calc_peakRH_linearize( - ds, - altitude="height", - rh="rh", - time_dim="sounding", - z_min=200.0, - z_max=900.0, - z_min_lin=50.0, - z_max_lin=400.0, -): - """ - Calculate maximum in relative humidity (RH) profile that minimizes difference between - observed RH profile and RH profile linearized around lower levels (i.e. 200-400m). - Assume moisture is well-mixed at these lower levels, such that RH would increase and - assume linear profile as an idealization - Inputs: - -- ds: dataset - -- altitude, rh, time_dim: variable names - -- z_min and z_max: lower and upper bounds for mixed layer height - -- z_min_lin and z_max_lin: bounds for linearization of observed RH profile - Outputs: - -- da: datarray containing h_peakRH_linfit - """ - h_peakRH_linfit = np.zeros(len(ds[rh])) - - dz = int(ds[altitude].diff(dim=altitude)[1]) # dz=10m - - mixed_layer = np.logical_and( - ds[altitude] >= z_min, ds[altitude] <= 1500 - ) # enforce z_max later - ml_linfit = np.logical_and( - ds[altitude] >= z_min_lin, ds[altitude] <= z_max_lin - ) # for linearized RH profile - - for i in range(len(ds[rh])): - - rh_profile = ds[rh].isel({time_dim: i}).interpolate_na(dim=altitude) - X_height = rh_profile[ml_linfit][altitude] - X = sm.add_constant(X_height.values) # add intercept - model = sm.OLS( - rh_profile[ml_linfit].values, X - ).fit() # instantiate linear model - linearized_RH = model.predict(sm.add_constant(ds[altitude])) - - # 'find_peaks' is scipy function - idx_peaks_RH_raw, _ = find_peaks(rh_profile[mixed_layer]) - # add +X index from mixed layer (ml) lower bound (lb), divided by dz=10m - add_to_index_ml_lb = int(z_min / dz) - idx_peaks_RH = idx_peaks_RH_raw + add_to_index_ml_lb - obs_RH_peaks = rh_profile[idx_peaks_RH] # RH values at the local maxima - linearized_RH_values = linearized_RH[ - idx_peaks_RH - ] # linearized RH values at the observed RH maxima - min_diff_RH = min( - abs(obs_RH_peaks - linearized_RH_values) - ) # find minimum RH difference between observed peaks and linearized profile - h_peakRH_linfit[i] = min_diff_RH[altitude] - - # set 'anomalously' high values > z_max to nan - # ... somewhat arbitrary choice - h_peakRH_linfit[h_peakRH_linfit > z_max] = np.nan - - dims = list(ds.dims.keys()) - # drop the height coord - del dims[dims.index(altitude)] - da = xr.DataArray(h_peakRH_linfit, dims=dims, coords={d: ds[d] for d in dims}) - da.attrs["long_name"] = "mixed layer height (from RH peak and linearization)" - da.attrs["units"] = "m" - return da - - -#%% callfunction - -z_peakRH_linfit = calc_peakRH_linearize( - all_sondes, - altitude="height", - rh="rh", - time_dim="sounding", - z_min=200.0, - z_max=900.0, - z_min_lin=50.0, - z_max_lin=200.0, -) - -#%% basic plots - -import seaborn as sns -import matplotlib.pyplot as plt -import numpy as np - -sns.set( - context="notebook", - style="whitegrid", - palette="deep", - font="sans-serif", - font_scale=2, - color_codes=True, - rc=None, -) -plt.figure(figsize=(15, 10)) -sns.distplot(z_peakRH_linfit) -plt.ylabel("pdf") -plt.xlabel("mixed layer height from RHmax linearize (m)") - -plt.figure(figsize=(15, 10)) -plt.scatter(np.arange(len(z_peakRH_linfit)), z_peakRH_linfit) -plt.xlabel("sonde number") -plt.ylabel("mixed layer height from RHmax linearize (m)") diff --git a/eurec4a_environment/variables/calc_density.py b/eurec4a_environment/variables/calc_density.py index bbda8ca..c2569fe 100644 --- a/eurec4a_environment/variables/calc_density.py +++ b/eurec4a_environment/variables/calc_density.py @@ -3,9 +3,8 @@ import numpy as np import xarray as xr - -# from ...constants import Rd, Rv, eps -from constants import Rd, Rv, eps +from ...constants import Rd, Rv, eps +#from constants import Rd, Rv, eps def calc_density(ds, pres="p", temp="T", specific_humidity="q", altitude="height"): From 1436af2947e6aafa5a06d1ecac65b1925d2ddc15 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Fri, 31 Jul 2020 16:44:09 +0200 Subject: [PATCH 12/17] Update eurec4a_environment/variables/calc_density.py Co-authored-by: Geet George <46085279+Geet-George@users.noreply.github.com> --- eurec4a_environment/variables/calc_density.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/eurec4a_environment/variables/calc_density.py b/eurec4a_environment/variables/calc_density.py index c2569fe..353a819 100644 --- a/eurec4a_environment/variables/calc_density.py +++ b/eurec4a_environment/variables/calc_density.py @@ -3,7 +3,7 @@ import numpy as np import xarray as xr -from ...constants import Rd, Rv, eps +from ..constants import Rd, Rv, eps #from constants import Rd, Rv, eps From 5924a80353e2fb4737956f9c65a226585c6db188 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Wed, 12 Aug 2020 09:31:13 +0200 Subject: [PATCH 13/17] Delete mixed_layer_height_from_gradient.py superfluous file --- .../mixed_layer_height_from_gradient.py | 88 ------------------- 1 file changed, 88 deletions(-) delete mode 100644 eurec4a_environment/variables/boundary_layer/mixed_layer_height_from_gradient.py diff --git a/eurec4a_environment/variables/boundary_layer/mixed_layer_height_from_gradient.py b/eurec4a_environment/variables/boundary_layer/mixed_layer_height_from_gradient.py deleted file mode 100644 index 49868a9..0000000 --- a/eurec4a_environment/variables/boundary_layer/mixed_layer_height_from_gradient.py +++ /dev/null @@ -1,88 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Find mixed layer height as layer over which x(z+1) - x(z) < threshold - -Inputs: - -- ds: Dataset - -- var: variable name - -- threshold: in units of variable (i.e. g/kg or K). Maximum difference allowed in one vertical step, dz=10m, i.e. threshold_q = 0.4 g/kg, - or vertical gradient dq/dz < 0.04 g/kg/m - -Outputs: DataArray containing mixed layer height from gradient method - -Note that function is quite slow and should be optimized - -""" - -import numpy as np -import xarray as xr -from scipy.signal import find_peaks -import statsmodels.api as sm -from constants import Rd, Rv, eps -import seaborn as sns - -# import eurec4a_environment.source_data -# ds = eurec4a_environment.source_data.open_joanne_dataset() - - -def calc_from_gradient( - ds, var, threshold, z_min=200, altitude="height", time="sounding" -): - def calculateHmix_var( - density_profile, - var_profile, - threshold, - z_min, - altitude="height", - time="sounding", - ): - - var_diff = 0 - numer = 0 - denom = 0 - var_mix = 0 - - k = int(z_min / 10) # enforce lower limit = 200m. k is index of lower limit - - # var(z) - weighted_mean(z-1) < threshold - while abs(var_diff) < threshold: - - numer += 0.5 * ( - density_profile[k + 1] * var_profile[k + 1] - + density_profile[k] * var_profile[k] - ) - denom += 0.5 * (density_profile[k + 1] + density_profile[k]) - var_mix = numer / denom - k += 1 - var_diff = (var_profile[k] - var_mix).values - - hmix = var_profile[altitude].values[k] - return hmix - - # call inner function - from eurec4a_environment.variables.calc_density import calc_density - - da_density = calc_density(ds) - - hmix_vec = np.zeros(len(ds[var])) - for i in range((len(ds[var]))): - density_profile = da_density.isel({time: i}) - var_profile = ds[var].isel({time: i}) - hmix_vec[i] = calculateHmix_var( - density_profile, - var_profile, - threshold, - z_min, - altitude="height", - time="sounding", - ) - - dims = list(ds.dims.keys()) - # drop the height coord - del dims[dims.index(altitude)] - da = xr.DataArray(hmix_vec, dims=dims, coords={d: ds[d] for d in dims}) - da.attrs["long_name"] = f"mixed layer height (from gradient method using {var})" - da.attrs["units"] = "m" - - return da From 13c4d4f5c56e83207d11d047b5e551430b63533c Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Thu, 20 Aug 2020 14:19:35 +0200 Subject: [PATCH 14/17] refactor remaining functions --- eurec4a_environment/plot/profile.py | 2 - .../boundary_layer/mixed_layer_height.py | 11 +- .../variables/calc_static_stability.py | 42 +++++++ .../inversion_layer_height.py | 116 ++++++++++++++++++ tests/test_boundary_layer.py | 18 ++- tests/test_inversion_layers.py | 25 ++++ tests/test_static_stability.py | 15 +++ 7 files changed, 218 insertions(+), 11 deletions(-) create mode 100755 eurec4a_environment/variables/calc_static_stability.py create mode 100755 eurec4a_environment/variables/inversion_layers/inversion_layer_height.py create mode 100755 tests/test_inversion_layers.py create mode 100755 tests/test_static_stability.py diff --git a/eurec4a_environment/plot/profile.py b/eurec4a_environment/plot/profile.py index d79c325..3dc252b 100644 --- a/eurec4a_environment/plot/profile.py +++ b/eurec4a_environment/plot/profile.py @@ -1,9 +1,7 @@ -import numpy as np import matplotlib.pyplot as plt import xarray as xr import matplotlib.dates as mdates import seaborn as sns -import typhon sns.set( context="notebook", diff --git a/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py b/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py index b4950c3..3dceda9 100644 --- a/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py +++ b/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py @@ -4,7 +4,7 @@ import statsmodels.api as sm -def calc_peak_RH(ds, altitude="alt", rh="RH", z_min=200.0, z_max=900.0): +def calc_peak_RH(ds, altitude="height", rh="rh", z_min=200.0, z_max=900.0): """ Calculate height at maximum relative humidity values """ @@ -28,7 +28,7 @@ def calc_peakRH_linearize( ds, altitude="height", rh="rh", - time_dim="sounding", + time="sounding", z_min=200.0, z_max=1500.0, z_min_lin=50.0, @@ -41,7 +41,7 @@ def calc_peakRH_linearize( assume linear profile as an idealization Inputs: -- ds: dataset - -- altitude, rh, time_dim: variable names + -- altitude, rh, time: variable names -- z_min and z_max: lower and upper bounds for mixed layer height -- z_min_lin and z_max_lin: bounds for linearization of observed RH profile Outputs: @@ -58,7 +58,7 @@ def calc_peakRH_linearize( for i in range(len(ds[rh])): - rh_profile = ds[rh].isel({time_dim: i}).interpolate_na(dim=altitude) + rh_profile = ds[rh].isel({time: i}).interpolate_na(dim=altitude) X_height = rh_profile[ml_linfit][altitude] X = sm.add_constant(X_height.values) # add intercept model = sm.OLS( @@ -107,8 +107,7 @@ def calc_from_gradient( Outputs: DataArray containing mixed layer height from gradient method - Note that function is quite slow and should be optimized - + Note that function is slow and should be optimized """ def calculateHmix_var( diff --git a/eurec4a_environment/variables/calc_static_stability.py b/eurec4a_environment/variables/calc_static_stability.py new file mode 100755 index 0000000..74257bf --- /dev/null +++ b/eurec4a_environment/variables/calc_static_stability.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import xarray as xr + + +def calc_static_stability(ds, pres="p", temp="T", theta="theta", altitude="height"): + + """ Equation from: https://www.ncl.ucar.edu/Document/Functions/Contributed/static_stability.shtml + S = -(T/theta)*d(theta)/dp. From Bluestein (1992), pg 197, eqn 4.3.8 + + Static stability measures the gravitational resistance of an atmosphere to vertical displacements. + It results from fundamental buoyant adjustments, and so it is determined by the vertical stratification + of density or potential temperature. + """ + + # !! edit these lines based on units function !! + # Celsius to Kelvin + if ds[temp].max().values < 100: + temp_K = ds[temp] + 273.15 + else: + temp_K = ds[temp] + + # convert Pa to hPa + if ds[pres].max().values > 1400: + ds[pres] = ds[pres] / 100 + + static_stability = -(temp_K / ds[theta]) * ( + ds[theta].diff(dim=altitude) / ds[pres].diff(dim=altitude) + ) + static_stability = static_stability.transpose( + transpose_coords=True + ) # same dimension order as original dataset + + dims = list(static_stability.dims) + da = xr.DataArray( + static_stability, dims=dims, coords={d: static_stability[d] for d in dims} + ) + da.attrs["long_name"] = "static stability" + da.attrs["units"] = "K/hPa" + + return da diff --git a/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py b/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py new file mode 100755 index 0000000..cfb45e4 --- /dev/null +++ b/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np +import xarray as xr + + +def calc_inversion_static_stability( + ds, time="sounding", altitude="height", z_min=1000, z_max=3000 +): + + """ + Returns inversion height defined as the height where static stability profiles + first exceed 0.1 K/hPa between given heights + Adapted from Bony and Stevens 2018, Measuring Area-Averaged Vertical Motions with Dropsondes, p 772 + To do: combine static stability criterion with one for moisture (i.e. hydrolapse) + """ + + from eurec4a_environment.variables.calc_static_stability import ( + calc_static_stability, + ) + + ds_inversion = ds.sel({altitude: slice(z_min, z_max)}) + da_static_stability = calc_static_stability(ds_inversion) + + inversion_height_stability = np.zeros(len(da_static_stability[time])) + for i in range(len(da_static_stability[time])): + stability_sounding = da_static_stability.isel({time: i}).values + idx_stability = np.argmax(stability_sounding > 0.1) + inversion_height_stability[i] = da_static_stability[altitude][idx_stability] + + dims = list(da_static_stability.dims) + # drop the height coord + del dims[dims.index(altitude)] + + da = xr.DataArray( + inversion_height_stability, + dims=dims, + coords={d: da_static_stability[d] for d in dims}, + ) + da.attrs["long_name"] = "inversion base height, static stability > 0.1 K/hPa" + da.attrs["units"] = "m" + + return da + + +def calc_inversion_grad_RH_T( + ds, + time="sounding", + altitude="height", + temperature="T", + rh="rh", + z_min=1000, + z_max=2500, +): + + """ + Inversion base following Grindinger et al, 1992; Cao et al 2007 + dT/dz > 0 and dRH/dz < 0 + To do: find all inversion bases and heights, not just first inversion base + """ + + ds_inversion = ds.sel({altitude: slice(z_min, z_max)}) + gradient_RH = ds_inversion[rh].differentiate(coord=altitude) + gradient_T = ds_inversion[temperature].differentiate(coord=altitude) + inversion_levels = gradient_T[altitude].where( + (gradient_RH <= 0) & (gradient_T >= 0) + ) + + num_iters = len(inversion_levels[time]) + inv_base_height_gradients = np.zeros(num_iters) + for i in range(num_iters): + alt_inv_da = inversion_levels.isel({time: i}) + inversion_alts = alt_inv_da[~np.isnan(alt_inv_da)] + if inversion_alts[altitude].size == 0: + print(f"sounding {i} has no inversion below {z_max} m") + continue + inv_base_height_gradients[i] = inversion_alts[0] # find first inversion base + inv_base_height_gradients[inv_base_height_gradients == 0] = np.nan + + dims = list(ds.dims.keys()) + # drop the height coord + del dims[dims.index(altitude)] + da = xr.DataArray( + inv_base_height_gradients, dims=dims, coords={d: ds[d] for d in dims} + ) + da.attrs["long_name"] = "inversion base height from RH and T gradient" + da.attrs["units"] = "m" + return da + + # sns.distplot(inv_base_height_gradients) + + # plot vertical profiles dRH/dz and dT/dz + + # fig, ax = plt.subplots(1,2,figsize=(20,8)) + # for i in range(5): # len(gradient_RH) + # ax[0].plot(gradient_RH.isel({time: i}), gradient_RH[altitude].values,color="lightgrey", linewidth=2,alpha=0.5) + # ax[0].plot(gradient_RH.mean(dim=time),gradient_RH[altitude].values,linewidth=4, color='black') + # ax[0].spines['right'].set_visible(False) + # ax[0].spines['top'].set_visible(False) + # ax[0].set_xlabel('$\partial$RH/$\partial$z') + # ax[0].set_ylabel('Altitude / m') + # ax[0].set_ylim([0,z_max]) + # ax[0].axvline(x=0, color='black') + # ax[0].set_xlim([ -1, 1]) + + # for i in range(5): + # ax[1].plot(gradient_T.isel({time: i}), gradient_T[altitude].values,color="lightgrey", linewidth=2,alpha=0.5) + # ax[1].plot(gradient_T.mean(dim=time),gradient_T[altitude].values,linewidth=4, color='black') + # ax[1].spines['right'].set_visible(False) + # ax[1].spines['top'].set_visible(False) + # ax[1].set_xlabel('$\partial$T/$\partial$z (dz=10m)' ) + # #ax[1].set_ylabel('Altitude / m') + # ax[1].set_ylim([0,z_max]) + # ax[1].axvline(x=0, color='black') + # ax[1].set_xlim([-0.03, 0.03]) diff --git a/tests/test_boundary_layer.py b/tests/test_boundary_layer.py index a1acaa8..85a3f42 100644 --- a/tests/test_boundary_layer.py +++ b/tests/test_boundary_layer.py @@ -15,9 +15,9 @@ def test_LCL_Bolton(ds_isentropic_test_profiles): def test_mixed_layer_height_RHmax(ds_isentropic_test_profiles): ds = ds_isentropic_test_profiles # set the RH profile so we know where the peak should be - z0 = 600. + z0 = 600.0 z = ds.alt - ds['rh'] = 1.0 - np.maximum(np.abs(z0 - z), 0) / z0 + ds["rh"] = 1.0 - np.maximum(np.abs(z0 - z), 0) / z0 da_rh_peak = boundary_layer.mixed_layer_height.calc_peak_RH( ds=ds, altitude="alt", rh="rh" ) @@ -29,4 +29,16 @@ def test_mixed_layer_height_RH_lin(): ds = ds.isel(sounding=slice(0, 10)) da_lin = boundary_layer.mixed_layer_height.calc_peakRH_linearize(ds) assert len(da_lin) == 10 - assert np.all(da_lin < 1500.) + assert np.all(da_lin < 1500.0) + + +def test_mixed_layer_height_gradient(): + ds = eurec4a_environment.source_data.open_joanne_dataset() + ds = ds.isel(sounding=slice(0, 10)) + # adjust threshold based on units 0.4 g/kg/10m or 0.0004 kg/kg/10m + da_gradient = boundary_layer.mixed_layer_height.calc_from_gradient( + ds, var="q", threshold=0.4 / 1000 + ) + assert len(da_gradient) == 10 + assert np.all(da_gradient < 1500.0) + assert da_gradient.units == "m" diff --git a/tests/test_inversion_layers.py b/tests/test_inversion_layers.py new file mode 100755 index 0000000..23ff466 --- /dev/null +++ b/tests/test_inversion_layers.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np +from eurec4a_environment.variables import inversion_layers +import eurec4a_environment.source_data + + +def test_inversion_static_stability(): + ds = eurec4a_environment.source_data.open_joanne_dataset() + ds = ds.isel(sounding=slice(0, 10)) + da_ss = inversion_layers.inversion_layer_height.calc_inversion_static_stability(ds) + assert len(da_ss) == 10 + assert np.all(da_ss < 5000.0) + assert da_ss.units == "m" + + +def test_inversion_gradient_RH_T(): + ds = eurec4a_environment.source_data.open_joanne_dataset() + ds = ds.isel(sounding=slice(0, 10)) + # adjust threshold based on units 0.4 g/kg/10m or 0.0004 kg/kg/10m + da_gradRH_T = inversion_layers.inversion_layer_height.calc_inversion_grad_RH_T(ds) + assert len(da_gradRH_T) == 10 + assert np.all(da_gradRH_T < 5000.0) + assert da_gradRH_T.units == "m" diff --git a/tests/test_static_stability.py b/tests/test_static_stability.py new file mode 100755 index 0000000..2a6e6f1 --- /dev/null +++ b/tests/test_static_stability.py @@ -0,0 +1,15 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + + +from eurec4a_environment.variables.calc_static_stability import calc_static_stability +import eurec4a_environment.source_data + + +def test_static_stability_calculation(): + ds = eurec4a_environment.source_data.open_joanne_dataset() + da_static_stability = calc_static_stability(ds) + # check that we get some sensible numbers + assert da_static_stability.units == "K/hPa" + assert da_static_stability.max() < 3 + assert da_static_stability.min() > -3 From c0bbb7c4be4b1f0cfe300dc3d0e9c34372509ca6 Mon Sep 17 00:00:00 2001 From: Anna Lea Albright Date: Mon, 24 Aug 2020 11:46:36 +0200 Subject: [PATCH 15/17] minor changes --- .../variables/boundary_layer/lcl.py | 2 +- .../inversion_layer_height.py | 27 +------------------ tests/test_boundary_layer.py | 2 +- tests/test_inversion_layers.py | 6 ++--- 4 files changed, 6 insertions(+), 31 deletions(-) diff --git a/eurec4a_environment/variables/boundary_layer/lcl.py b/eurec4a_environment/variables/boundary_layer/lcl.py index 5d2f836..c081784 100644 --- a/eurec4a_environment/variables/boundary_layer/lcl.py +++ b/eurec4a_environment/variables/boundary_layer/lcl.py @@ -3,7 +3,7 @@ from ...constants import cp_d, g -def find_LCL_Bolton(ds, temperature="T", rh="RH", altitude="alt"): +def find_LCL_Bolton(ds, temperature="T", rh="rh", altitude="height"): """ Calculates distribution of LCL from RH and T at different vertical levels returns mean LCL from this distribution diff --git a/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py b/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py index cfb45e4..e28db0c 100755 --- a/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py +++ b/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py @@ -88,29 +88,4 @@ def calc_inversion_grad_RH_T( da.attrs["units"] = "m" return da - # sns.distplot(inv_base_height_gradients) - - # plot vertical profiles dRH/dz and dT/dz - - # fig, ax = plt.subplots(1,2,figsize=(20,8)) - # for i in range(5): # len(gradient_RH) - # ax[0].plot(gradient_RH.isel({time: i}), gradient_RH[altitude].values,color="lightgrey", linewidth=2,alpha=0.5) - # ax[0].plot(gradient_RH.mean(dim=time),gradient_RH[altitude].values,linewidth=4, color='black') - # ax[0].spines['right'].set_visible(False) - # ax[0].spines['top'].set_visible(False) - # ax[0].set_xlabel('$\partial$RH/$\partial$z') - # ax[0].set_ylabel('Altitude / m') - # ax[0].set_ylim([0,z_max]) - # ax[0].axvline(x=0, color='black') - # ax[0].set_xlim([ -1, 1]) - - # for i in range(5): - # ax[1].plot(gradient_T.isel({time: i}), gradient_T[altitude].values,color="lightgrey", linewidth=2,alpha=0.5) - # ax[1].plot(gradient_T.mean(dim=time),gradient_T[altitude].values,linewidth=4, color='black') - # ax[1].spines['right'].set_visible(False) - # ax[1].spines['top'].set_visible(False) - # ax[1].set_xlabel('$\partial$T/$\partial$z (dz=10m)' ) - # #ax[1].set_ylabel('Altitude / m') - # ax[1].set_ylim([0,z_max]) - # ax[1].axvline(x=0, color='black') - # ax[1].set_xlim([-0.03, 0.03]) + diff --git a/tests/test_boundary_layer.py b/tests/test_boundary_layer.py index 85a3f42..45a2c49 100644 --- a/tests/test_boundary_layer.py +++ b/tests/test_boundary_layer.py @@ -19,7 +19,7 @@ def test_mixed_layer_height_RHmax(ds_isentropic_test_profiles): z = ds.alt ds["rh"] = 1.0 - np.maximum(np.abs(z0 - z), 0) / z0 da_rh_peak = boundary_layer.mixed_layer_height.calc_peak_RH( - ds=ds, altitude="alt", rh="rh" + ds=ds, altitude="height", rh="rh" ) assert np.allclose(da_rh_peak, z0) diff --git a/tests/test_inversion_layers.py b/tests/test_inversion_layers.py index 23ff466..5620210 100755 --- a/tests/test_inversion_layers.py +++ b/tests/test_inversion_layers.py @@ -2,14 +2,14 @@ # -*- coding: utf-8 -*- import numpy as np -from eurec4a_environment.variables import inversion_layers +from eurec4a_environment.variables.inversion_layers import inversion_layer_height import eurec4a_environment.source_data def test_inversion_static_stability(): ds = eurec4a_environment.source_data.open_joanne_dataset() ds = ds.isel(sounding=slice(0, 10)) - da_ss = inversion_layers.inversion_layer_height.calc_inversion_static_stability(ds) + da_ss = inversion_layer_height.calc_inversion_static_stability(ds) assert len(da_ss) == 10 assert np.all(da_ss < 5000.0) assert da_ss.units == "m" @@ -19,7 +19,7 @@ def test_inversion_gradient_RH_T(): ds = eurec4a_environment.source_data.open_joanne_dataset() ds = ds.isel(sounding=slice(0, 10)) # adjust threshold based on units 0.4 g/kg/10m or 0.0004 kg/kg/10m - da_gradRH_T = inversion_layers.inversion_layer_height.calc_inversion_grad_RH_T(ds) + da_gradRH_T = inversion_layer_height.calc_inversion_grad_RH_T(ds) assert len(da_gradRH_T) == 10 assert np.all(da_gradRH_T < 5000.0) assert da_gradRH_T.units == "m" From b619ec1fc35218339bd2d3dd9d53aa96c62dc06b Mon Sep 17 00:00:00 2001 From: Leif Denby Date: Wed, 3 Feb 2021 20:31:56 +0000 Subject: [PATCH 16/17] apply black and fix flake8 errors --- eurec4a_environment/__init__.py | 4 +++- eurec4a_environment/plot/profile.py | 7 ++++--- eurec4a_environment/unit.py | 2 -- .../variables/boundary_layer/__init__.py | 2 +- .../variables/boundary_layer/inversion_height.py | 2 -- .../variables/boundary_layer/lcl.py | 5 ++++- .../variables/calc_static_stability.py | 15 +++++++++------ .../inversion_layers/inversion_layer_height.py | 2 -- eurec4a_environment/variables/utils.py | 2 +- tests/test_boundary_layer.py | 5 +++-- tests/test_tropical.py | 4 +++- 11 files changed, 28 insertions(+), 22 deletions(-) diff --git a/eurec4a_environment/__init__.py b/eurec4a_environment/__init__.py index c2d5f4d..52f9b57 100644 --- a/eurec4a_environment/__init__.py +++ b/eurec4a_environment/__init__.py @@ -37,7 +37,9 @@ def get_fields(ds, *names, **names_and_units): dataarrays = [] coords = {} - all_fields = itertools.chain(names_and_units.items(), zip(names, [None] * len(names))) + all_fields = itertools.chain( + names_and_units.items(), zip(names, [None] * len(names)) + ) for (name, units) in all_fields: da = get_field(ds=ds, name=name, units=units) diff --git a/eurec4a_environment/plot/profile.py b/eurec4a_environment/plot/profile.py index 99c41a7..baa3b70 100644 --- a/eurec4a_environment/plot/profile.py +++ b/eurec4a_environment/plot/profile.py @@ -1,5 +1,4 @@ import matplotlib.pyplot as plt -import xarray as xr import matplotlib.dates as mdates import seaborn as sns @@ -91,8 +90,10 @@ def plot_profile_1D( height_labels=None, **kwargs ): - """Plot a 1D profile plot of two variables and a height level - input a 1D dataset (that is profiles e.g. of one timestep or sounding or a mean profile) + """ + Plot a 1D profile plot of two variables and a height level input a 1D + dataset (that is profiles e.g. of one timestep or sounding or a mean + profile) """ fig, ax1 = plt.subplots(figsize=(12, 12)) diff --git a/eurec4a_environment/unit.py b/eurec4a_environment/unit.py index 9afff24..dcee7e8 100644 --- a/eurec4a_environment/unit.py +++ b/eurec4a_environment/unit.py @@ -5,8 +5,6 @@ HAS_UDUNITS2 = True except FileNotFoundError: - import warnings - HAS_UDUNITS2 = False diff --git a/eurec4a_environment/variables/boundary_layer/__init__.py b/eurec4a_environment/variables/boundary_layer/__init__.py index 5ed4712..f5ea189 100644 --- a/eurec4a_environment/variables/boundary_layer/__init__.py +++ b/eurec4a_environment/variables/boundary_layer/__init__.py @@ -1,2 +1,2 @@ from . import lcl # noqa -from . import mixed_layer_height +from . import mixed_layer_height # noqa diff --git a/eurec4a_environment/variables/boundary_layer/inversion_height.py b/eurec4a_environment/variables/boundary_layer/inversion_height.py index 2d87863..cafa0f5 100644 --- a/eurec4a_environment/variables/boundary_layer/inversion_height.py +++ b/eurec4a_environment/variables/boundary_layer/inversion_height.py @@ -1,5 +1,3 @@ -import xarray as xr - from ... import get_field from ... import nomenclature as nom diff --git a/eurec4a_environment/variables/boundary_layer/lcl.py b/eurec4a_environment/variables/boundary_layer/lcl.py index ec10d6b..2405b97 100644 --- a/eurec4a_environment/variables/boundary_layer/lcl.py +++ b/eurec4a_environment/variables/boundary_layer/lcl.py @@ -7,7 +7,10 @@ def find_LCL_Bolton( - ds, temperature=nom.TEMPERATURE, rh=nom.RELATIVE_HUMIDITY, altitude=nom.ALTITUDE, + ds, + temperature=nom.TEMPERATURE, + rh=nom.RELATIVE_HUMIDITY, + altitude=nom.ALTITUDE, vertical_coord=nom.ALTITUDE, ): """ diff --git a/eurec4a_environment/variables/calc_static_stability.py b/eurec4a_environment/variables/calc_static_stability.py index 74257bf..e62002e 100755 --- a/eurec4a_environment/variables/calc_static_stability.py +++ b/eurec4a_environment/variables/calc_static_stability.py @@ -5,13 +5,16 @@ def calc_static_stability(ds, pres="p", temp="T", theta="theta", altitude="height"): + """ + Equation from: + https://www.ncl.ucar.edu/Document/Functions/Contributed/static_stability.shtml + + S = -(T/theta)*d(theta)/dp. From Bluestein (1992), pg 197, eqn 4.3.8 - """ Equation from: https://www.ncl.ucar.edu/Document/Functions/Contributed/static_stability.shtml - S = -(T/theta)*d(theta)/dp. From Bluestein (1992), pg 197, eqn 4.3.8 - - Static stability measures the gravitational resistance of an atmosphere to vertical displacements. - It results from fundamental buoyant adjustments, and so it is determined by the vertical stratification - of density or potential temperature. + Static stability measures the gravitational resistance of an atmosphere to + vertical displacements. It results from fundamental buoyant adjustments, + and so it is determined by the vertical stratification of density or + potential temperature. """ # !! edit these lines based on units function !! diff --git a/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py b/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py index e28db0c..154e390 100755 --- a/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py +++ b/eurec4a_environment/variables/inversion_layers/inversion_layer_height.py @@ -87,5 +87,3 @@ def calc_inversion_grad_RH_T( da.attrs["long_name"] = "inversion base height from RH and T gradient" da.attrs["units"] = "m" return da - - diff --git a/eurec4a_environment/variables/utils.py b/eurec4a_environment/variables/utils.py index f37ed4c..1d06528 100644 --- a/eurec4a_environment/variables/utils.py +++ b/eurec4a_environment/variables/utils.py @@ -1,5 +1,5 @@ def apply_by_column(ds, vertical_coord, fn): - if not vertical_coord in ds.dims: + if vertical_coord not in ds.dims: raise Exception( f"`{vertical_coord}` is not a dimension of the provided " "dataset" ) diff --git a/tests/test_boundary_layer.py b/tests/test_boundary_layer.py index f328b07..a88b239 100644 --- a/tests/test_boundary_layer.py +++ b/tests/test_boundary_layer.py @@ -3,6 +3,7 @@ from eurec4a_environment.variables import boundary_layer from eurec4a_environment.variables.boundary_layer import inversion_height from eurec4a_environment import nomenclature as nom +import eurec4a_environment.source_data def test_LCL_Bolton(ds_isentropic_test_profiles): @@ -44,7 +45,7 @@ def test_mixed_layer_height_gradient(): assert np.all(da_gradient < 1500.0) assert da_gradient.units == "m" - + def test_inversion_height_gradient_RH(ds_isentropic_test_profiles): ds = ds_isentropic_test_profiles.copy() z_INV = 2000.0 @@ -54,4 +55,4 @@ def test_inversion_height_gradient_RH(ds_isentropic_test_profiles): da_inv = boundary_layer.inversion_height.find_inversion_height_grad_RH( ds=ds, rh=nom.RELATIVE_HUMIDITY ) - assert np.allclose(da_inv, z_INV, atol=20) ## within 20m \ No newline at end of file + assert np.allclose(da_inv, z_INV, atol=20) ## within 20m diff --git a/tests/test_tropical.py b/tests/test_tropical.py index ed070a7..d4966d2 100644 --- a/tests/test_tropical.py +++ b/tests/test_tropical.py @@ -99,7 +99,9 @@ def test_eis_joanne_profile(): ) ds[LTS_name] = da_lts - da_eis = tropical_variables.estimated_inversion_strength(ds=ds, LCL=LCL_name, LTS=LTS_name) + da_eis = tropical_variables.estimated_inversion_strength( + ds=ds, LCL=LCL_name, LTS=LTS_name + ) assert da_eis.units == "K" assert np.all(da_eis > 0.0) From cbd6bb379c7d48ca573821be94f01f8f1ca405c1 Mon Sep 17 00:00:00 2001 From: Leif Denby Date: Wed, 3 Feb 2021 21:12:36 +0000 Subject: [PATCH 17/17] start refactoring mixed-layer peak RH linfit calculation --- .../boundary_layer/mixed_layer_height.py | 118 ++++++++++++------ 1 file changed, 78 insertions(+), 40 deletions(-) diff --git a/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py b/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py index 62b4fd0..aa6e16e 100644 --- a/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py +++ b/eurec4a_environment/variables/boundary_layer/mixed_layer_height.py @@ -3,7 +3,10 @@ from scipy.signal import find_peaks import statsmodels.api as sm -from ... import get_field, nomenclature as nom +from ... import get_field, get_fields, nomenclature as nom +from ..utils import apply_by_column + +DEBUG = False def calc_peak_RH( @@ -37,36 +40,45 @@ def calc_peakRH_linearize( z_max=1500.0, z_min_lin=50.0, z_max_lin=400.0, + vertical_coord=":altitude", ): """ - Calculate maximum in relative humidity (RH) profile that minimizes difference between - observed RH profile and RH profile linearized around lower levels (i.e. 200-400m). - Assume moisture is well-mixed at these lower levels, such that RH would increase and - assume linear profile as an idealization - Inputs: - -- ds: dataset - -- altitude, rh, time: variable names - -- z_min and z_max: lower and upper bounds for mixed layer height - -- z_min_lin and z_max_lin: bounds for linearization of observed RH profile - Outputs: - -- da: datarray containing h_peakRH_linfit - """ - da_rh = get_field(ds=ds, name=rh, units="%") - da_alt = get_field(ds=ds, name=altitude, units="m") - - h_peakRH_linfit = np.zeros(len(da_rh)) - - dz = int(da_alt.diff(dim=altitude)[1]) # dz=10m - - mixed_layer = np.logical_and(da_alt >= z_min, da_alt <= z_max) - - ml_linfit = np.logical_and( - da_alt >= z_min_lin, da_alt <= z_max_lin - ) # for linearized RH profile - - for i in range(len(ds[rh])): - - rh_profile = ds[rh].isel({time: i}).interpolate_na(dim=altitude) + Calculate maximum in relative humidity (RH) profile that minimizes + difference between observed RH profile and RH profile linearized around + lower levels (i.e. 200-400m). Assume moisture is well-mixed at these lower + levels, such that RH would increase and assume linear profile as an + idealization. + + Inputs: + - ds: dataset + - altitude, rh, time: variable names + - z_min and z_max: lower and upper bounds for mixed layer height + - z_min_lin and z_max_lin: bounds for linearization of observed RH profile + + Outputs: + - da: datarray containing h_peakRH_linfit + """ + # if `vertical_coord` hasn't been changed we'll assume that the vertical + # coordinate has name of the altitude variable + if vertical_coord == ":altitude": + vertical_coord = altitude + + import matplotlib.pyplot as plt + + def _calc_mixed_layer_depth(ds_column): + da_alt = ds_column[altitude] + + mixed_layer = np.logical_and(da_alt >= z_min, da_alt <= z_max) + ml_linfit = np.logical_and( + da_alt >= z_min_lin, da_alt <= z_max_lin + ) # for linearized RH profile + + dz = int(da_alt.diff(dim=altitude)[1]) # dz=10m + + rh_profile = ds_column[rh].interpolate_na(dim=altitude) + + # fit a straight line to the relative humidity values in the + # "ml-linfit" region X_height = rh_profile[ml_linfit][altitude] X = sm.add_constant(X_height.values) # add intercept model = sm.OLS( @@ -79,6 +91,7 @@ def calc_peakRH_linearize( # add +X index from mixed layer (ml) lower bound (lb), divided by dz=10m add_to_index_ml_lb = int(z_min / dz) idx_peaks_RH = idx_peaks_RH_raw + add_to_index_ml_lb + obs_RH_peaks = rh_profile[idx_peaks_RH] # RH values at the local maxima linearized_RH_values = linearized_RH[ idx_peaks_RH @@ -86,19 +99,44 @@ def calc_peakRH_linearize( min_diff_RH = min( abs(obs_RH_peaks - linearized_RH_values) ) # find minimum RH difference between observed peaks and linearized profile - h_peakRH_linfit[i] = min_diff_RH[altitude] - # set 'anomalously' high values > z_max to nan - # ... somewhat arbitrary choice - h_peakRH_linfit[h_peakRH_linfit > z_max] = np.nan + h_peakRH_linfit = min_diff_RH[altitude] - dims = list(ds.dims.keys()) - # drop the height coord - del dims[dims.index(altitude)] - da = xr.DataArray(h_peakRH_linfit, dims=dims, coords={d: ds[d] for d in dims}) - da.attrs["long_name"] = "mixed layer height (from RH peak and linearization)" - da.attrs["units"] = "m" - return da + if DEBUG: + ds_column[rh].plot(y=altitude, label="raw profile") + rh_profile.plot(y=altitude, label="raw profile interpolated") + plt.plot(linearized_RH, da_alt, label="RH interpolation in linfit region") + + rh_profile.isel({altitude: idx_peaks_RH}).plot( + marker="x", linestyle="--", label="interpolated profile peaks", y=altitude + ) + + plt.axhline(h_peakRH_linfit, label="h_peakRH_linfit") + + plt.legend() + plt.show() + + # set 'anomalously' high values > z_max to nan + # ... somewhat arbitrary choice + if h_peakRH_linfit > z_max: + h_peakRH_linfit = np.nan + + import ipdb + + ipdb.set_trace() + + return h_peakRH_linfit + + ds_subset = get_fields(ds=ds, **{rh: "%", altitude: "m"}) + da_ml_depth = apply_by_column( + ds=ds_subset, vertical_coord=vertical_coord, fn=_calc_mixed_layer_depth + ) + + da_ml_depth.attrs[ + "long_name" + ] = "mixed layer height (from RH peak and linearization)" + da_ml_depth.attrs["units"] = "m" + return da_ml_depth def calc_from_gradient(