Remove cloud cover, SHF/LHF and upper limit on surface temperature fo…

…r bedrock stations

#155

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']) 
 

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

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'):