Smoothing and inter/extrapolating gps coordinates

.github/workflows/process_test.yml

@@ -38,7 +38,7 @@ jobs:
         run: |
           mkdir $GITHUB_WORKSPACE/out/
           for i in $(echo ${{ env.TEST_STATION }} | tr ' ' '\n'); do
-            python3 $GITHUB_WORKSPACE/main/src/pypromice/process/get_l3.py -v $GITHUB_WORKSPACE/main/src/pypromice/process/variables.csv -m $GITHUB_WORKSPACE/main/src/pypromice/process/metadata.csv -c $GITHUB_WORKSPACE/aws-l0/raw/config/$i.toml -i $GITHUB_WORKSPACE/aws-l0/raw -o $GITHUB_WORKSPACE/out/
+            python3 $GITHUB_WORKSPACE/main/src/pypromice/process/get_l2.py -v $GITHUB_WORKSPACE/main/src/pypromice/process/variables.csv -m $GITHUB_WORKSPACE/main/src/pypromice/process/metadata.csv -c $GITHUB_WORKSPACE/aws-l0/raw/config/$i.toml -i $GITHUB_WORKSPACE/aws-l0/raw -o $GITHUB_WORKSPACE/out/
           done
       - name: Upload test output
         uses: actions/upload-artifact@v3

.gitignore

@@ -29,3 +29,4 @@ src/pypromice/postprocess/positions.csv
 
 # sqlite db files
 *.db
+*.bak

setup.py

@@ -42,8 +42,9 @@
     'console_scripts': [
         'get_promice_data = pypromice.get.get_promice_data:get_promice_data',
         'get_l0tx = pypromice.tx.get_l0tx:get_l0tx',
+        'get_l2 = pypromice.process.get_l2:get_l2',
+        'join_l2 = pypromice.process.join_l2:join_l2',
         'get_l3 = pypromice.process.get_l3:get_l3',
-        'join_l3 = pypromice.process.join_l3:join_l3',
         'get_watsontx = pypromice.tx.get_watsontx:get_watsontx',
         'get_bufr = pypromice.postprocess.get_bufr:main',
         'get_msg = pypromice.tx.get_msg:get_msg'

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

@@ -30,7 +30,18 @@ def toL2(
     eps_clear=9.36508e-6,
     emissivity=0.97,
 ) -> xr.Dataset:
-    '''Process one Level 1 (L1) product to Level 2
+    '''Process one Level 1 (L1) product to Level 2.
+    In this step we do:
+        - manual flagging and adjustments
+        - automated QC: persistence, percentile
+        - 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
+        - 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
+        - calculation of directional wind speed
 
     Parameters
     ----------
@@ -59,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
@@ -85,10 +97,20 @@ def toL2(
     ds['dlr'] = ds.dlr.where(ds.t_rad.notnull())
     ds['ulr'] = ds.ulr.where(ds.t_rad.notnull())
 
+    # calculating realtive humidity with regard to ice
     T_100 = _getTempK(T_0)
     ds['rh_u_cor'] = correctHumidity(ds['rh_u'], ds['t_u'],
                                      T_0, T_100, ews, ei0)
 
+    if ds.attrs['number_of_booms']==2:
+        ds['rh_l_cor'] = correctHumidity(ds['rh_l'], ds['t_l'],
+                                         T_0, T_100, ews, ei0)
+
+    if hasattr(ds,'t_i'):
+        if ~ds['t_i'].isnull().all():
+            ds['rh_i_cor'] = correctHumidity(ds['rh_i'], ds['t_i'],
+                                             T_0, T_100, ews, ei0)
+
     # Determiune cloud cover for on-ice stations
     cc = calcCloudCoverage(ds['t_u'], T_0, eps_overcast, eps_clear,        # Calculate cloud coverage
                            ds['dlr'], ds.attrs['station_id'])
@@ -176,22 +198,52 @@ def toL2(
         ds['precip_u_cor'], ds['precip_u_rate'] = correctPrecip(ds['precip_u'],
                                                                 ds['wspd_u'])
     if ds.attrs['number_of_booms']==2:
-        ds['rh_l_cor'] = correctHumidity(ds['rh_l'], ds['t_l'],           # Correct relative humidity
-                                         T_0, T_100, ews, ei0)
-
         if ~ds['precip_l'].isnull().all() and precip_flag:                     # Correct precipitation
             ds['precip_l_cor'], ds['precip_l_rate']= correctPrecip(ds['precip_l'],
                                                                    ds['wspd_l'])
 
-    if hasattr(ds,'t_i'):
-        if ~ds['t_i'].isnull().all():                                          # Instantaneous msg processing
-            ds['rh_i_cor'] = correctHumidity(ds['rh_i'], ds['t_i'],       # Correct relative humidity
-                                             T_0, T_100, ews, ei0)
+    # Get directional wind speed 
+    ds['wdir_u'] = ds['wdir_u'].where(ds['wspd_u'] != 0)                   
+    ds['wspd_x_u'], ds['wspd_y_u'] = calcDirWindSpeeds(ds['wspd_u'], ds['wdir_u']) 
+
+    if ds.attrs['number_of_booms']==2:                                         
+        ds['wdir_l'] = ds['wdir_l'].where(ds['wspd_l'] != 0) 
+        ds['wspd_x_l'], ds['wspd_y_l'] = calcDirWindSpeeds(ds['wspd_l'], ds['wdir_l'])
+
+    if hasattr(ds, 'wdir_i'):
+        if ~ds['wdir_i'].isnull().all() and ~ds['wspd_i'].isnull().all():
+            ds['wdir_i'] = ds['wdir_i'].where(ds['wspd_i'] != 0)                   
+            ds['wspd_x_i'], ds['wspd_y_i'] = calcDirWindSpeeds(ds['wspd_i'], ds['wdir_i'])   
+
 
     ds = clip_values(ds, vars_df)
     return ds
 
 
+def calcDirWindSpeeds(wspd, wdir, deg2rad=np.pi/180):
+    '''Calculate directional wind speed from wind speed and direction
+
+    Parameters
+    ----------
+    wspd : xr.Dataarray
+        Wind speed data array
+    wdir : xr.Dataarray
+        Wind direction data array
+    deg2rad : float
+        Degree to radians coefficient. The default is np.pi/180
+
+    Returns
+    -------
+    wspd_x : xr.Dataarray
+        Wind speed in X direction
+    wspd_y : xr.Datarray
+        Wind speed in Y direction
+    '''        
+    wspd_x = wspd * np.sin(wdir * deg2rad)
+    wspd_y = wspd * np.cos(wdir * deg2rad) 
+    return wspd_x, wspd_y
+
+
 def calcCloudCoverage(T, T_0, eps_overcast, eps_clear, dlr, station_id):
     '''Calculate cloud cover from T and T_0