brought changes from #252

src/pypromice/process/L0toL1.py

@@ -5,9 +5,9 @@
 import numpy as np
 import pandas as pd
 import xarray as xr
-import re
-
+import re, logging
 from pypromice.process.value_clipping import clip_values
+logger = logging.getLogger(__name__)
 
 
 def toL1(L0, vars_df, T_0=273.15, tilt_threshold=-100):
@@ -28,9 +28,10 @@ def toL1(L0, vars_df, T_0=273.15, tilt_threshold=-100):
     -------
     ds : xarray.Dataset
         Level 1 dataset
-    '''
+    '''    
     assert(type(L0) == xr.Dataset)
     ds = L0
+    ds.attrs['level'] = 'L1'
 
     for l in list(ds.keys()):
         if l not in ['time', 'msg_i', 'gps_lat', 'gps_lon', 'gps_alt', 'gps_time']:
@@ -247,7 +248,7 @@ def getPressDepth(z_pt, p, pt_antifreeze, pt_z_factor, pt_z_coef, pt_z_p_coef):
         rho_af = 1145
     else:
         rho_af = np.nan
-        print('ERROR: Incorrect metadata: "pt_antifreeze" = ' +
+        logger.info('ERROR: Incorrect metadata: "pt_antifreeze" = ' +
               f'{pt_antifreeze}. Antifreeze mix only supported at 50% or 100%')
         # assert(False)
 

src/pypromice/process/L1toL2.py

@@ -37,11 +37,11 @@ def toL2(
         - custom filter: gps_alt filter, NaN t_rad removed from dlr & ulr
         - smoothing of tilt and rot
         - calculation of rh with regards to ice in subfreezin conditions
-        - caluclation of cloud coverage
+        - calculation of cloud coverage
         - correction of dsr and usr for tilt
         - filtering of dsr based on a theoritical TOA irradiance and grazing light
         - calculation of albedo
-        - caluclation of directional wind speed
+        - calculation of directional wind speed
 
     Parameters
     ----------
@@ -70,6 +70,7 @@ def toL2(
         Level 2 dataset
     '''
     ds = L1.copy(deep=True)                                                    # Reassign dataset
+    ds.attrs['level'] = 'L2'
     try:
         ds = adjustTime(ds)                                                    # Adjust time after a user-defined csv files
         ds = flagNAN(ds)                                                       # Flag NaNs after a user-defined csv files

src/pypromice/process/L2toL3.py

@@ -23,6 +23,7 @@ def toL3(L2, T_0=273.15):
         Steam point temperature. Default is 273.15.
     '''
     ds = L2
+    ds.attrs['level'] = 'L3'
 
     T_100 = _getTempK(T_0)                                                     # Get steam point temperature as K 
 
@@ -123,75 +124,6 @@ def toL3(L2, T_0=273.15):
     return ds
 
 
-def find_breaks(df,alpha):
-    '''Detects potential relocation of the station from the GPS measurements.
-    The code first makes a forward linear interpolation of the coordinates and
-    then looks for important jumps in latitude, longitude and altitude. The jumps
-    that are higher than a given threshold (expressed as a multiple of the 
-    standard deviation) are mostly caused by the station being moved during
-    maintenance. To avoid misclassification, only the jumps detected in May-Sept. 
-    are kept.
-    
-    Parameters
-    ----------
-    df : pandas.Series
-        series of observed latitude, longitude or elevation
-    alpha: float
-        coefficient to be applied to the the standard deviation of the daily
-        coordinate fluctuation
-    '''
-    diff = df.resample('D').median().interpolate(
-        method='linear', limit_area='inside', limit_direction='forward').diff()        
-    thresh = diff.std() * alpha
-    list_diff = diff.loc[diff.abs()>thresh].reset_index()
-    list_diff = list_diff.loc[list_diff.time.dt.month.isin([5,6,7,8,9])]
-    list_diff['year']=list_diff.time.dt.year
-    list_diff=list_diff.groupby('year').max()
-    return diff, [None]+list_diff.time.to_list()+[None]
-
-
-def piecewise_smoothing_and_interpolation(df_in, breaks):
-    '''Smoothes, inter- or extrapolate the gps observations. The processing is 
-    done piecewise so that each period between station relocation are done 
-    separately (no smoothing of the jump due to relocation). Locally Weighted
-    Scatterplot Smoothing (lowess) is then used to smooth the available
-    observations. Then this smoothed curve is interpolated linearly over internal
-    gaps. Eventually, this interpolated curve is extrapolated linearly for 
-    timestamps before the first valid measurement and after the last valid
-    measurement.
-    
-    Parameters
-    ----------
-    df_in : pandas.Series
-        series of observed latitude, longitude or elevation
-    breaks: list
-        List of timestamps of station relocation. First and last item should be
-        None so that they can be used in slice(breaks[i], breaks[i+1])
-    '''
-    df_all = pd.Series() # dataframe gathering all the smoothed pieces
-    for i in range(len(breaks)-1):
-        df = df_in.loc[slice(breaks[i], breaks[i+1])].copy()
-
-        y_sm = lowess(df,
-                      pd.to_numeric(df.index),
-                      is_sorted=True, frac=1/3, it=0,
-                      )
-        df.loc[df.notnull()] = y_sm[:,1]
-        df = df.interpolate(method='linear', limit_area='inside')
-
-        last_valid_6_months = slice(df.last_valid_index()-pd.to_timedelta('180D'),None)
-        df.loc[last_valid_6_months] = (df.loc[last_valid_6_months].interpolate( axis=0,
-            method='spline',order=1, limit_direction='forward', fill_value="extrapolate")).values
-
-        first_valid_6_months = slice(None, df.first_valid_index()+pd.to_timedelta('180D'))
-        df.loc[first_valid_6_months] = (df.loc[first_valid_6_months].interpolate( axis=0,
-            method='spline',order=1, limit_direction='backward', fill_value="extrapolate")).values
-        df_all=pd.concat((df_all, df))
-
-    df_all = df_all[~df_all.index.duplicated(keep='first')]
-    return df_all.values  
-
-
 def calcHeatFlux(T_0, T_h, Tsurf_h, WS_h, z_WS, z_T, q_h, p_h, 
                 kappa=0.4, WS_lim=1., z_0=0.001, g=9.82, es_0=6.1071, eps=0.622, 
                 gamma=16., L_sub=2.83e6, L_dif_max=0.01, c_pd=1005., aa=0.7, 
@@ -439,6 +371,74 @@ def calcSpHumid(T_0, T_100, T_h, p_h, RH_cor_h, es_0=6.1071, es_100=1013.246, ep
     return RH_cor_h * q_sat / 100 
 
 
+def find_breaks(df,alpha):
+    '''Detects potential relocation of the station from the GPS measurements.
+    The code first makes a forward linear interpolation of the coordinates and
+    then looks for important jumps in latitude, longitude and altitude. The jumps
+    that are higher than a given threshold (expressed as a multiple of the 
+    standard deviation) are mostly caused by the station being moved during
+    maintenance. To avoid misclassification, only the jumps detected in May-Sept. 
+    are kept.
+    
+    Parameters
+    ----------
+    df : pandas.Series
+        series of observed latitude, longitude or elevation
+    alpha: float
+        coefficient to be applied to the the standard deviation of the daily
+        coordinate fluctuation
+    '''
+    diff = df.resample('D').median().interpolate(
+        method='linear', limit_area='inside', limit_direction='forward').diff()        
+    thresh = diff.std() * alpha
+    list_diff = diff.loc[diff.abs()>thresh].reset_index()
+    list_diff = list_diff.loc[list_diff.time.dt.month.isin([5,6,7,8,9])]
+    list_diff['year']=list_diff.time.dt.year
+    list_diff=list_diff.groupby('year').max()
+    return diff, [None]+list_diff.time.to_list()+[None]
+
+
+def piecewise_smoothing_and_interpolation(df_in, breaks):
+    '''Smoothes, inter- or extrapolate the gps observations. The processing is 
+    done piecewise so that each period between station relocation are done 
+    separately (no smoothing of the jump due to relocation). Locally Weighted
+    Scatterplot Smoothing (lowess) is then used to smooth the available
+    observations. Then this smoothed curve is interpolated linearly over internal
+    gaps. Eventually, this interpolated curve is extrapolated linearly for 
+    timestamps before the first valid measurement and after the last valid
+    measurement.
+    
+    Parameters
+    ----------
+    df_in : pandas.Series
+        series of observed latitude, longitude or elevation
+    breaks: list
+        List of timestamps of station relocation. First and last item should be
+        None so that they can be used in slice(breaks[i], breaks[i+1])
+    '''
+    df_all = pd.Series() # dataframe gathering all the smoothed pieces
+    for i in range(len(breaks)-1):
+        df = df_in.loc[slice(breaks[i], breaks[i+1])].copy()
+
+        y_sm = lowess(df,
+                      pd.to_numeric(df.index),
+                      is_sorted=True, frac=1/3, it=0,
+                      )
+        df.loc[df.notnull()] = y_sm[:,1]
+        df = df.interpolate(method='linear', limit_area='inside')
+
+        last_valid_6_months = slice(df.last_valid_index()-pd.to_timedelta('180D'),None)
+        df.loc[last_valid_6_months] = (df.loc[last_valid_6_months].interpolate( axis=0,
+            method='spline',order=1, limit_direction='forward', fill_value="extrapolate")).values
+
+        first_valid_6_months = slice(None, df.first_valid_index()+pd.to_timedelta('180D'))
+        df.loc[first_valid_6_months] = (df.loc[first_valid_6_months].interpolate( axis=0,
+            method='spline',order=1, limit_direction='backward', fill_value="extrapolate")).values
+        df_all=pd.concat((df_all, df))
+
+    df_all = df_all[~df_all.index.duplicated(keep='first')]
+    return df_all.values
+
 def _calcAtmosDens(p_h, R_d, T_h, T_0):                                        # TODO: check this shouldn't be in this step somewhere
     '''Calculate atmospheric density'''
     return 100 * p_h / R_d / (T_h + T_0)