From 7e615347fa2cbe02f78359805183937b12220331 Mon Sep 17 00:00:00 2001 From: BaptisteVandecrux Date: Thu, 10 Aug 2023 11:10:30 +0200 Subject: [PATCH] Remove cloud cover, SHF/LHF and upper limit on surface temperature for bedrock stations #155 --- src/pypromice/process/L0toL1.py | 11 +++++++---- src/pypromice/process/L1toL2.py | 15 +++++++++------ src/pypromice/process/L2toL3.py | 30 ++++++++++++++++-------------- 3 files changed, 32 insertions(+), 24 deletions(-) diff --git a/src/pypromice/process/L0toL1.py b/src/pypromice/process/L0toL1.py index 8e953642..b1d66d30 100644 --- a/src/pypromice/process/L0toL1.py +++ b/src/pypromice/process/L0toL1.py @@ -8,7 +8,7 @@ import re -def toL1(L0, vars_df, flag_file=None, T_0=273.15, tilt_threshold=-100): +def toL1(L0, vars_df, T_0=273.15, tilt_threshold=-100): '''Process one Level 0 (L0) product to Level 1 Parameters @@ -17,8 +17,6 @@ def toL1(L0, vars_df, flag_file=None, T_0=273.15, tilt_threshold=-100): Level 0 dataset vars_df : pd.DataFrame Metadata dataframe - flag_file : str - Flag .csv file path for bad data T_0 : int Air temperature for sonic ranger adjustment tilt_threshold : int @@ -91,9 +89,14 @@ def toL1(L0, vars_df, flag_file=None, T_0=273.15, tilt_threshold=-100): if hasattr(ds, 'bedrock'): # Fix tilt to zero if station is on bedrock if ds.attrs['bedrock']==True or ds.attrs['bedrock'].lower() in 'true': + ds.attrs['bedrock'] = True # ensures all AWS objects have a 'bedrock' attribute ds['tilt_x'] = (('time'), np.arange(ds['time'].size)*0) ds['tilt_y'] = (('time'), np.arange(ds['time'].size)*0) - + else: + ds.attrs['bedrock'] = False # ensures all AWS objects have a 'bedrock' attribute + else: + ds.attrs['bedrock'] = False # ensures all AWS objects have a 'bedrock' attribute + ds['wdir_u'] = ds['wdir_u'].where(ds['wspd_u'] != 0) # Get directional wind speed ds['wspd_x_u'], ds['wspd_y_u'] = calcWindDir(ds['wspd_u'], ds['wdir_u']) diff --git a/src/pypromice/process/L1toL2.py b/src/pypromice/process/L1toL2.py index 625f1408..e5f9c63a 100644 --- a/src/pypromice/process/L1toL2.py +++ b/src/pypromice/process/L1toL2.py @@ -50,15 +50,18 @@ def toL2(L1, T_0=273.15, ews=1013.246, ei0=6.1071, eps_overcast=1., ds['rh_u_cor'] = correctHumidity(ds['rh_u'], ds['t_u'], T_0, T_100, ews, ei0) - # Determiune cloud cover - cc = calcCloudCoverage(ds['t_u'], T_0, eps_overcast, eps_clear, # Calculate cloud coverage - ds['dlr'], ds.attrs['station_id']) - ds['cc'] = (('time'), cc.data) + # Determiune cloud cover for on-ice stations + if not ds.attrs['bedrock']: + cc = calcCloudCoverage(ds['t_u'], T_0, eps_overcast, eps_clear, # Calculate cloud coverage + ds['dlr'], ds.attrs['station_id']) + ds['cc'] = (('time'), cc.data) # Determine surface temperature ds['t_surf'] = calcSurfaceTemperature(T_0, ds['ulr'], ds['dlr'], # Calculate surface temperature emissivity) - + if not ds.attrs['bedrock']: + ds['t_surf'] = ds['t_surf'].where(ds['t_surf'] <= 0, other = 0) + # Determine station position relative to sun doy = ds['time'].to_dataframe().index.dayofyear.values # Gather variables to calculate sun pos hour = ds['time'].to_dataframe().index.hour.values @@ -489,7 +492,7 @@ def calcSurfaceTemperature(T_0, ulr, dlr, emissivity): Calculated surface temperature ''' t_surf = ((ulr - (1 - emissivity) * dlr) / emissivity / 5.67e-8)**0.25 - T_0 - return t_surf.where(t_surf <= 0, other = 0) + return t_surf def calcTilt(tilt_x, tilt_y, deg2rad): '''Calculate station tilt diff --git a/src/pypromice/process/L2toL3.py b/src/pypromice/process/L2toL3.py index d005cb42..66f7b345 100755 --- a/src/pypromice/process/L2toL3.py +++ b/src/pypromice/process/L2toL3.py @@ -45,15 +45,16 @@ def toL3(L2, T_0=273.15, z_0=0.001, R_d=287.05, eps=0.622, es_0=6.1071, rho_atm_u = 100 * p_h_u / R_d / (T_h_u + T_0) # Calculate atmospheric density nu_u = calcVisc(T_h_u, T_0, rho_atm_u) # Calculate kinematic viscosity q_h_u = calcHumid(T_0, T_100, T_h_u, es_0, es_100, eps, # Calculate specific humidity - p_h_u, RH_cor_h_u) - SHF_h_u, LHF_h_u= calcHeatFlux(T_0, T_h_u, Tsurf_h, rho_atm_u, WS_h_u, # Calculate latent and sensible heat fluxes - z_WS_u, z_T_u, nu_u, q_h_u, p_h_u) - SHF_h_u, LHF_h_u = cleanHeatFlux(SHF_h_u, LHF_h_u, T_h_u, Tsurf_h, p_h_u, # Clean heat flux values - WS_h_u, RH_cor_h_u, ds['z_boom_u']) + p_h_u, RH_cor_h_u) + if not ds.attrs['bedrock']: + SHF_h_u, LHF_h_u= calcHeatFlux(T_0, T_h_u, Tsurf_h, rho_atm_u, WS_h_u, # Calculate latent and sensible heat fluxes + z_WS_u, z_T_u, nu_u, q_h_u, p_h_u) + SHF_h_u, LHF_h_u = cleanHeatFlux(SHF_h_u, LHF_h_u, T_h_u, Tsurf_h, p_h_u, # Clean heat flux values + WS_h_u, RH_cor_h_u, ds['z_boom_u']) + ds['dshf_u'] = (('time'), SHF_h_u.data) + ds['dlhf_u'] = (('time'), LHF_h_u.data) q_h_u = 1000 * q_h_u # Convert sp.humid from kg/kg to g/kg q_h_u = cleanSpHumid(q_h_u, T_h_u, Tsurf_h, p_h_u, RH_cor_h_u) # Clean sp.humid values - ds['dshf_u'] = (('time'), SHF_h_u.data) - ds['dlhf_u'] = (('time'), LHF_h_u.data) ds['qh_u'] = (('time'), q_h_u.data) # Lower boom bulk calculation @@ -70,15 +71,16 @@ def toL3(L2, T_0=273.15, z_0=0.001, R_d=287.05, eps=0.622, es_0=6.1071, rho_atm_l = 100 * p_h_l / R_d / (T_h_l + T_0) # Calculate atmospheric density nu_l = calcVisc(T_h_l, T_0, rho_atm_l) # Calculate kinematic viscosity q_h_l = calcHumid(T_0, T_100, T_h_l, es_0, es_100, eps, # Calculate sp.humidity - p_h_l, RH_cor_h_l) - SHF_h_l, LHF_h_l= calcHeatFlux(T_0, T_h_l, Tsurf_h, rho_atm_l, WS_h_l, # Calculate latent and sensible heat fluxes - z_WS_l, z_T_l, nu_l, q_h_l, p_h_l) - SHF_h_l, LHF_h_l = cleanHeatFlux(SHF_h_l, LHF_h_l, T_h_l, Tsurf_h, p_h_l, # Clean heat flux values - WS_h_l, RH_cor_h_l, ds['z_boom_l']) + p_h_l, RH_cor_h_l) + if not ds.attrs['bedrock']: + SHF_h_l, LHF_h_l= calcHeatFlux(T_0, T_h_l, Tsurf_h, rho_atm_l, WS_h_l, # Calculate latent and sensible heat fluxes + z_WS_l, z_T_l, nu_l, q_h_l, p_h_l) + SHF_h_l, LHF_h_l = cleanHeatFlux(SHF_h_l, LHF_h_l, T_h_l, Tsurf_h, p_h_l, # Clean heat flux values + WS_h_l, RH_cor_h_l, ds['z_boom_l']) + ds['dshf_l'] = (('time'), SHF_h_l.data) + ds['dlhf_l'] = (('time'), LHF_h_l.data) q_h_l = 1000 * q_h_l # Convert sp.humid from kg/kg to g/kg q_h_l = cleanSpHumid(q_h_l, T_h_l, Tsurf_h, p_h_l, RH_cor_h_l) # Clean sp.humid values - ds['dshf_l'] = (('time'), SHF_h_l.data) - ds['dlhf_l'] = (('time'), LHF_h_l.data) ds['qh_l'] = (('time'), q_h_l.data) # if hasattr(ds, 'wspd_x_i'):