Write L2 to file and leave all additional variable derivation for the L2toL3 step

.github/workflows/process_l2_test.yml

@@ -0,0 +1,52 @@
+on:
+  pull_request:
+    types: [opened, reopened, synchronize, edited]
+  workflow_dispatch:
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    name: process_l2_test
+    steps:
+      - name: Install Python
+        uses: actions/setup-python@v4
+        with:
+          python-version: "3.8"        
+      - name: Checkout repo
+        uses: actions/checkout@v3
+        with:
+          path: "main"
+          token: ${{ secrets.GITHUB_TOKEN }}
+      - name: Install dependencies
+        shell: bash
+        run: |
+          python -m pip install --upgrade pip
+          pip install wheel
+          python3 -m pip install --upgrade setuptools
+          cd $GITHUB_WORKSPACE/main
+          pip install .  
+      - name: Clone AWS Level 0 data repo for testing
+        env:
+          GITLAB_TOKEN : ${{ secrets.GITLAB_TOKEN }}
+        run: | 
+          cd $GITHUB_WORKSPACE
+          git clone --depth 1 https://oauth2:${{ env.GITLAB_TOKEN }}@geusgitlab.geus.dk/glaciology-and-climate/promice/aws-l0.git
+      - name: Run L0 to L3 processing
+        env:
+          TEST_STATION: KAN_U HUM
+        shell: bash
+        run: |
+          mkdir $GITHUB_WORKSPACE/out/
+          mkdir $GITHUB_WORKSPACE/out/L2/
+          for i in $(echo ${{ env.TEST_STATION }} | tr ' ' '\n'); do
+            python3 $GITHUB_WORKSPACE/main/src/pypromice/process/get_l2.py -c $GITHUB_WORKSPACE/aws-l0/tx/config/$i.toml -i $GITHUB_WORKSPACE/aws-l0/tx -o $GITHUB_WORKSPACE/out/L2/
+          done
+#          mkdir $GITHUB_WORKSPACE/out/L3/            
+#          for i in $(echo ${{ env.TEST_STATION }} | tr ' ' '\n'); do           
+#            python3 $GITHUB_WORKSPACE/main/src/pypromice/process/get_l2tol3.py -i $GITHUB_WORKSPACE/out/L2/$i/$i_hour.nc -o $GITHUB_WORKSPACE/out/ -t 60min
+#          done
+      - name: Upload test output
+        uses: actions/upload-artifact@v3
+        with:
+          name: result
+          path: out

.github/workflows/process_test.yml → .github/workflows/process_l3_test.yml

@@ -6,7 +6,7 @@ on:
 jobs:
   build:
     runs-on: ubuntu-latest
-    name: process_test
+    name: process_l3_test
     steps:
       - name: Install Python
         uses: actions/setup-python@v4
@@ -33,12 +33,12 @@ jobs:
           git clone --depth 1 https://oauth2:${{ env.GITLAB_TOKEN }}@geusgitlab.geus.dk/glaciology-and-climate/promice/aws-l0.git
       - name: Run data processing
         env:
-          TEST_STATION: KPC_U CEN2 JAR
+          TEST_STATION: KPC_U CEN2
         shell: bash
         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_l3.py -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

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="pypromice",
-    version="1.3.6",
+    version="1.4.0",
     author="GEUS Glaciology and Climate",
     description="PROMICE/GC-Net data processing toolbox",
     long_description=long_description,
@@ -36,14 +36,16 @@
         "pypromice.qc.percentiles": ["thresholds.csv"],
         "pypromice.postprocess": ["station_configurations.toml", "positions_seed.csv"],
     },
-    install_requires=['numpy>=1.23.0', 'pandas>=1.5.0', 'xarray>=2022.6.0', 'toml', 'scipy>=1.9.0', 'Bottleneck', 'netcdf4', 'pyDataverse', 'eccodes','scikit-learn>=1.1.0'],
+    install_requires=['numpy>=1.23.0', 'pandas>=1.5.0', 'xarray>=2022.6.0', 'toml', 'scipy>=1.9.0', 'Bottleneck', 'netcdf4', 'pyDataverse==0.3.1', 'eccodes', 'scikit-learn>=1.1.0'],
 #    extras_require={'postprocess': ['eccodes','scikit-learn>=1.1.0']},
     entry_points={
     'console_scripts': [
         'get_promice_data = pypromice.get.get_promice_data:get_promice_data',
         'get_l0tx = pypromice.tx.get_l0tx:get_l0tx',
+        'join_l2 = pypromice.process.join_l2:join_l2',
+        'get_l2 = pypromice.process.get_l2:get_l2',
         'get_l3 = pypromice.process.get_l3:get_l3',
-        'join_l3 = pypromice.process.join_l3:join_l3',
+        'get_l2tol3 = pypromice.process.get_l2tol3:get_l2tol3',
         '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/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
     ----------
@@ -85,10 +96,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 +197,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
 

src/pypromice/process/L2toL3.py

@@ -7,7 +7,10 @@
 
 def toL3(L2, T_0=273.15, z_0=0.001, R_d=287.05, eps=0.622, es_0=6.1071, 
          es_100=1013.246):
-    '''Process one Level 2 (L2) product to Level 3 (L3) 
+    '''Process one Level 2 (L2) product to Level 3 (L3) meaning calculating all
+    derived variables:
+        - Sensible fluxes
+
 
     Parameters
     ----------
@@ -32,9 +35,6 @@ def toL3(L2, T_0=273.15, z_0=0.001, R_d=287.05, eps=0.622, es_0=6.1071,
 
     T_100 = _getTempK(T_0)                                                     # Get steam point temperature as K 
 
-    ds['wdir_u'] = ds['wdir_u'].where(ds['wspd_u'] != 0)                       # Get directional wind speed                    
-    ds['wspd_x_u'], ds['wspd_y_u'] = calcDirWindSpeeds(ds['wspd_u'], ds['wdir_u']) 
-
     # Upper boom bulk calculation
     T_h_u = ds['t_u'].copy()                                                   # Copy for processing
     p_h_u = ds['p_u'].copy()
@@ -85,41 +85,9 @@ def toL3(L2, T_0=273.15, z_0=0.001, R_d=287.05, eps=0.622, es_0=6.1071,
         q_h_l = cleanSpHumid(q_h_l, T_h_l, Tsurf_h, p_h_l, RH_cor_h_l)         # Clean sp.humid values
         ds['qh_l'] = (('time'), q_h_l.data)    
 
-        ds['wdir_l'] = ds['wdir_l'].where(ds['wspd_l'] != 0)                   # Get directional wind speed    
-        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():      # Instantaneous msg processing
-            ds['wdir_i'] = ds['wdir_i'].where(ds['wspd_i'] != 0)               # Get directional wind speed                    
-            ds['wspd_x_i'], ds['wspd_y_i'] = calcDirWindSpeeds(ds['wspd_i'], ds['wdir_i'])   
-
     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 calcHeatFlux(T_0, T_h, Tsurf_h, rho_atm, WS_h, z_WS, z_T, nu, 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,