Fix OOL Dependency Handling, RH Clipping, and Bufr Export Rounding

src/pypromice/postprocess/bufr_utilities.py

@@ -120,7 +120,8 @@ class BUFRVariables:
     )
     # https://vocabulary-manager.eumetsat.int/vocabularies/BUFR/WMO/32/TABLE_B/012101
     # Scale: 2, unit: K
-    airTemperature: float = attrs.field(converter=round_converter(2))
+    # NOTE: The expected scale is 2, but our instantanous data is rounded to 1 decimal.
+    airTemperature: float = attrs.field(converter=round_converter(1))
     # There is also a Dewpoint temperature in this template: 012103 which is currently unused.
     # https://vocabulary-manager.eumetsat.int/vocabularies/BUFR/WMO/32/TABLE_B/012103
     # Scale: 0, unit: %

src/pypromice/process/L1toL2.py

@@ -94,7 +94,7 @@ def toL2(
                           .ffill().bfill())
     mask = (np.abs(ds.gps_alt - baseline_elevation) < 100) & ds.gps_alt.notnull()
     ds[['gps_alt','gps_lon', 'gps_lat']] = ds[['gps_alt','gps_lon', 'gps_lat']].where(mask)
-    
+
     # removing dlr and ulr that are missing t_rad
     # this is done now becasue t_rad can be filtered either manually or with persistence
     ds['dlr'] = ds.dlr.where(ds.t_rad.notnull())
@@ -113,12 +113,12 @@ def toL2(
         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'])
     ds['cc'] = (('time'), cc.data)
-    
+
     # Determine surface temperature
     ds['t_surf'] = calcSurfaceTemperature(T_0, ds['ulr'], ds['dlr'],           # Calculate surface temperature
                                           emissivity)
@@ -136,7 +136,7 @@ def toL2(
     else:
         lat = ds['gps_lat'].mean()
         lon = ds['gps_lon'].mean()
-    
+
     # smoothing tilt and rot
     ds['tilt_x'] = smoothTilt(ds['tilt_x'])
     ds['tilt_y'] = smoothTilt(ds['tilt_y'])
@@ -151,8 +151,8 @@ def toL2(
     ZenithAngle_rad, ZenithAngle_deg = calcZenith(lat, Declination_rad,        # Calculate zenith
                                                   HourAngle_rad, deg2rad,
                                                   rad2deg)
-    
-    
+
+
     # Correct Downwelling shortwave radiation
     DifFrac = 0.2 + 0.8 * cc
     CorFac_all = calcCorrectionFactor(Declination_rad, phi_sensor_rad,         # Calculate correction
@@ -186,9 +186,9 @@ def toL2(
     TOA_crit_nopass = (ds['dsr_cor'] > (0.9 * isr_toa + 10))                   # Determine filter
     ds['dsr_cor'][TOA_crit_nopass] = np.nan                                    # Apply filter and interpolate
     ds['usr_cor'][TOA_crit_nopass] = np.nan
-    
-    ds['dsr_cor'] = ds.dsr_cor.where(ds.dsr.notnull())  
-    ds['usr_cor'] = ds.usr_cor.where(ds.usr.notnull())  
+
+    ds['dsr_cor'] = ds.dsr_cor.where(ds.dsr.notnull())
+    ds['usr_cor'] = ds.usr_cor.where(ds.usr.notnull())
     # # Check sun position
     # sundown = ZenithAngle_deg >= 90
     # _checkSunPos(ds, OKalbedos, sundown, sunonlowerdome, TOA_crit_nopass)
@@ -205,27 +205,33 @@ def toL2(
             ds['precip_l_cor'], ds['precip_l_rate']= correctPrecip(ds['precip_l'],
                                                                    ds['wspd_l'])
 
-    # 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']) 
+    get_directional_wind_speed(ds)                                            # Get directional wind speed
+
+    ds = clip_values(ds, vars_df)
+    return ds
+
+def get_directional_wind_speed(ds: xr.Dataset) -> xr.Dataset:
+    """
+    Calculate directional wind speed from wind speed and direction and mutates the dataset
+    """
+
+    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) 
+    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)
+            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'])
     return ds
 
 
 def calcDirWindSpeeds(wspd, wdir, deg2rad=np.pi/180):
     '''Calculate directional wind speed from wind speed and direction
-    
+
     Parameters
     ----------
     wspd : xr.Dataarray
@@ -234,16 +240,16 @@ def calcDirWindSpeeds(wspd, wdir, deg2rad=np.pi/180):
         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) 
+    wspd_y = wspd * np.cos(wdir * deg2rad)
     return wspd_x, wspd_y
 
 
@@ -328,7 +334,7 @@ def smoothTilt(da: xr.DataArray, threshold=0.2):
         either X or Y smoothed tilt inclinometer measurements
     '''
     # we calculate the moving standard deviation over a 3-day sliding window
-    # hourly resampling is necessary to make sure the same threshold can be used 
+    # hourly resampling is necessary to make sure the same threshold can be used
     # for 10 min and hourly data
     moving_std_gap_filled = da.to_series().resample('h').median().rolling(
                     3*24, center=True, min_periods=2

src/pypromice/process/value_clipping.py

@@ -2,6 +2,8 @@
 import pandas
 import xarray
 
+from pypromice.utilities.dependency_graph import DependencyGraph
+
 
 def clip_values(
     ds: xarray.Dataset,
@@ -27,26 +29,37 @@ def clip_values(
     cols = ["lo", "hi", "OOL"]
     assert set(cols) <= set(var_configurations.columns)
 
-    variable_limits = var_configurations[cols].dropna(how="all")
-    for var, row in variable_limits.iterrows():
+    variable_limits = var_configurations[cols].assign(
+        dependents=lambda df: df.OOL.fillna("").str.split(),
+        # Find the closure of dependents using the DependencyGraph class
+        dependents_closure=lambda df: DependencyGraph.from_child_mapping(
+            df.dependents
+        ).child_closure_mapping(),
+    )
 
+    for var, row in variable_limits.iterrows():
         if var not in list(ds.variables):
             continue
 
-        if ~np.isnan(row.lo):
-            ds[var] = ds[var].where(ds[var] >= row.lo)
-        if ~np.isnan(row.hi):
-            ds[var] = ds[var].where(ds[var] <= row.hi)
-
-        other_vars = row.OOL
-        if isinstance(other_vars, str) and ~ds[var].isnull().all():
-            for o in other_vars.split():
-                if o not in list(ds.variables):
-                    continue
-                else:
-                    if ~np.isnan(row.lo):
-                        ds[var] = ds[var].where(ds[var] >= row.lo)
-                    if ~np.isnan(row.hi):
-                        ds[var] = ds[var].where(ds[var] <= row.hi)
+        # This is a special case for rh_u_cor and rh_l_cor where values are clipped to 0 and 100.
+        if var in ["rh_u_cor", "rh_l_cor"]:
+            # Nan inputs should stay nan
+            was_nan = ds[var].isnull()
+            if ~np.isnan(row.lo):
+                ds[var] = ds[var].where(ds[var] >= row.lo, other=0)
+            if ~np.isnan(row.hi):
+                ds[var] = ds[var].where( ds[var] <= row.hi, other=100)
+            ds[var] = ds[var].where(~was_nan)
+        else:
+            if ~np.isnan(row.lo):
+                ds[var] = ds[var].where(ds[var] >= row.lo)
+            if ~np.isnan(row.hi):
+                ds[var] = ds[var].where(ds[var] <= row.hi)
+
+        # Flag dependents as NaN if parent is NaN
+        for o in row.dependents_closure:
+            if o not in list(ds.variables):
+                continue
+            ds[o] = ds[o].where(ds[var].notnull())
 
     return ds

src/pypromice/resources/variables.csv

@@ -11,15 +11,15 @@ qh_u,specific_humidity,Specific humidity (upper boom),kg/kg,modelResult,time,FAL
 rh_l,relative_humidity,Relative humidity (lower boom),%,physicalMeasurement,time,FALSE,,0,100,rh_l_cor,two-boom,1,1,1,4
 rh_l_cor,relative_humidity_corrected,Relative humidity (lower boom) - corrected,%,modelResult,time,FALSE,L2 or later,0,150,dshf_l dlhf_l qh_l,two-boom,0,1,1,4
 qh_l,specific_humidity,Specific humidity (lower boom),kg/kg,modelResult,time,FALSE,L2 or later,0,100,,two-boom,0,1,1,4
-wspd_u,wind_speed,Wind speed (upper boom),m s-1,physicalMeasurement,time,FALSE,,0,100,"wdir_u wspd_x_u wspd_y_u dshf_u dlhf_u qh_u, precip_u",all,1,1,1,4
-wspd_l,wind_speed,Wind speed (lower boom),m s-1,physicalMeasurement,time,FALSE,,0,100,"wdir_l wspd_x_l wspd_y_l dshf_l dlhf_l qh_l , precip_l",two-boom,1,1,1,4
+wspd_u,wind_speed,Wind speed (upper boom),m s-1,physicalMeasurement,time,FALSE,,0,100,wdir_u wspd_x_u wspd_y_u dshf_u dlhf_u qh_u precip_u,all,1,1,1,4
+wspd_l,wind_speed,Wind speed (lower boom),m s-1,physicalMeasurement,time,FALSE,,0,100,wdir_l wspd_x_l wspd_y_l dshf_l dlhf_l qh_l precip_l,two-boom,1,1,1,4
 wdir_u,wind_from_direction,Wind from direction (upper boom),degrees,physicalMeasurement,time,FALSE,,1,360,wspd_x_u wspd_y_u,all,1,1,1,4
 wdir_std_u,wind_from_direction_standard_deviation,Wind from direction (standard deviation),degrees,qualityInformation,time,FALSE,L0 or L2,,,,one-boom,1,1,0,4
 wdir_l,wind_from_direction,Wind from direction (lower boom),degrees,physicalMeasurement,time,FALSE,,1,360,wspd_x_l wspd_y_l,two-boom,1,1,1,4
-wspd_x_u,wind_speed_from_x_direction,Wind speed from x direction (upper boom),m s-1,modelResult,time,FALSE,L0 or L2,-100,100,wdir_u wspd_u,all,0,1,1,4
-wspd_y_u,wind_speed_from_y_direction,Wind speed from y direction (upper boom),m s-1,modelResult,time,FALSE,L0 or L2,-100,100,wdir_u wspd_u,all,0,1,1,4
-wspd_x_l,wind_speed_from_x_direction,Wind speed from x direction (lower boom),m s-1,modelResult,time,FALSE,L0 or L2,-100,100,wdir_l wspd_l,two-boom,0,1,1,4
-wspd_y_l,wind_speed_from_y_direction,Wind speed from y direction (lower boom),m s-1,modelResult,time,FALSE,L0 or L2,-100,100,wdir_l wspd_l,two-boom,0,1,1,4
+wspd_x_u,wind_speed_from_x_direction,Wind speed from x direction (upper boom),m s-1,modelResult,time,FALSE,L0 or L2,-100,100,"",all,0,1,1,4
+wspd_y_u,wind_speed_from_y_direction,Wind speed from y direction (upper boom),m s-1,modelResult,time,FALSE,L0 or L2,-100,100,"",all,0,1,1,4
+wspd_x_l,wind_speed_from_x_direction,Wind speed from x direction (lower boom),m s-1,modelResult,time,FALSE,L0 or L2,-100,100,"",two-boom,0,1,1,4
+wspd_y_l,wind_speed_from_y_direction,Wind speed from y direction (lower boom),m s-1,modelResult,time,FALSE,L0 or L2,-100,100,"",two-boom,0,1,1,4
 dsr,surface_downwelling_shortwave_flux,Downwelling shortwave radiation,W m-2,physicalMeasurement,time,FALSE,,-10,1500,albedo dsr_cor usr_cor,all,1,1,1,4
 dsr_cor,surface_downwelling_shortwave_flux_corrected,Downwelling shortwave radiation - corrected,W m-2,modelResult,time,FALSE,L2 or later,,,,all,0,1,1,4
 usr,surface_upwelling_shortwave_flux,Upwelling shortwave radiation,W m-2,physicalMeasurement,time,FALSE,,-10,1000,albedo dsr_cor usr_cor,all,1,1,1,4
@@ -102,5 +102,5 @@ rh_i,relative_humidity,Relative humidity (instantaneous),%,physicalMeasurement,t
 rh_i_cor,relative_humidity_corrected,Relative humidity (instantaneous) – corrected,%,modelResult,time,TRUE,L2 or later,0,100,,all,0,1,1,4
 wspd_i,wind_speed,Wind speed (instantaneous),m s-1,physicalMeasurement,time,TRUE,,0,100,wdir_i wspd_x_i wspd_y_i,all,1,1,1,4
 wdir_i,wind_from_direction,Wind from direction (instantaneous),degrees,physicalMeasurement,time,TRUE,,1,360,wspd_x_i wspd_y_i,all,1,1,1,4
-wspd_x_i,wind_speed_from_x_direction,Wind speed from x direction (instantaneous),m s-1,modelResult,time,TRUE,L2 or later,-100,100,wdir_i wspd_i,all,0,1,1,4
-wspd_y_i,wind_speed_from_y_direction,Wind speed from y direction (instantaneous),m s-1,modelResult,time,TRUE,L2 or later,-100,100,wdir_i wspd_i,all,0,1,1,4
+wspd_x_i,wind_speed_from_x_direction,Wind speed from x direction (instantaneous),m s-1,modelResult,time,TRUE,L2 or later,-100,100,"",all,0,1,1,4
+wspd_y_i,wind_speed_from_y_direction,Wind speed from y direction (instantaneous),m s-1,modelResult,time,TRUE,L2 or later,-100,100,"",all,0,1,1,4

src/pypromice/utilities/dependency_graph.py

@@ -0,0 +1,101 @@
+from typing import Mapping, Set, MutableMapping, Optional
+
+import attr
+
+__all__ = [
+    "DependencyNode",
+    "DependencyGraph",
+]
+
+
+@attr.define
+class DependencyNode:
+    name: str = attr.field()
+    parents: Set["DependencyNode"] = attr.field(factory=set)
+    children: Set["DependencyNode"] = attr.field(factory=set)
+
+    def add_child(self, child: "DependencyNode"):
+        self.children.add(child)
+        child.parents.add(self)
+
+    def get_children_closure(self, closure: Optional[Set[str]] = None) -> Set[str]:
+        is_root = closure is None
+        if closure is None:
+            closure = set()
+        if self.name in closure:
+            return closure
+        closure.add(self.name)
+        for child in self.children:
+            closure |= child.get_children_closure(closure)
+
+        if is_root:
+            closure.remove(self.name)
+        return closure
+
+    def get_parents_closure(self, closure: Optional[Set[str]] = None) -> Set[str]:
+        is_root = closure is None
+        if closure is None:
+            closure = set()
+        if self.name in closure:
+            return closure
+        closure.add(self.name)
+        for parent in self.parents:
+            closure |= parent.get_parents_closure(closure)
+
+        if is_root:
+            closure.remove(self.name)
+        return closure
+
+    def __hash__(self):
+        return hash(self.name)
+
+
+@attr.define
+class DependencyGraph:
+    nodes: MutableMapping[str, DependencyNode] = attr.field(factory=dict)
+
+    def add_node(self, name: str) -> DependencyNode:
+        if name not in self.nodes:
+            self.nodes[name] = DependencyNode(name)
+        return self.nodes[name]
+
+    def add_edge(self, parent: str, child: str):
+        parent_node = self.add_node(parent)
+        child_node = self.add_node(child)
+        parent_node.add_child(child_node)
+
+    @classmethod
+    def from_child_mapping(cls, mapping: Mapping[str, Set[str]]) -> "DependencyGraph":
+        graph = cls()
+        for var, children in mapping.items():
+            graph.add_node(var)
+            for child in children:
+                graph.add_edge(var, child)
+        return graph
+
+    @classmethod
+    def from_parent_mapping(cls, mapping: Mapping[str, Set[str]]) -> "DependencyGraph":
+        graph = cls()
+        for var, parents in mapping.items():
+            graph.add_node(var)
+            for parent in parents:
+                graph.add_edge(parent, var)
+        return graph
+
+    def child_mapping(self) -> Mapping[str, Set[str]]:
+        return {
+            node.name: {child.name for child in node.children}
+            for node in self.nodes.values()
+        }
+
+    def child_closure_mapping(self) -> Mapping[str, Set[str]]:
+        return {node.name: node.get_children_closure() for node in self.nodes.values()}
+
+    def parent_mapping(self) -> Mapping[str, Set[str]]:
+        return {
+            node.name: {parent.name for parent in node.parents}
+            for node in self.nodes.values()
+        }
+
+    def parent_closure_mapping(self) -> Mapping[str, Set[str]]:
+        return {node.name: node.get_parents_closure() for node in self.nodes.values()}