Merge pull request #156 from hgloeckner/bad_gpsalt_recovery

Bad gpsalt recovery

pydropsonde/circles.py

@@ -116,7 +116,7 @@ def get_xy_coords_for_circles(self):
         circle_altitude_attrs = {
             "long_name": "circle_altitude",
             "description": "Mean altitude of the aircraft during the circle",
-            "units": self.circle_ds.alt.attrs["units"],
+            "units": self.circle_ds[self.alt_dim].attrs["units"],
         }
         circle_time_attrs = {
             "long_name": "circle_time",

pydropsonde/helper/quality.py

@@ -99,10 +99,10 @@ def alt_below_aircraft(
         """
         alt_dim = self.alt_dim
         ds = self.qc_ds
-        self.qc_flags[f"{alt_dim}_below_aircraft"] = (
+        self.qc_flags["altitude_below_aircraft"] = (
             np.nanmax(ds[alt_dim].values) < maxalt
         )
-        if not self.qc_flags[f"{alt_dim}_below_aircraft"]:
+        if not self.qc_flags["altitude_below_aircraft"]:
             variables = ["lat", "lon", "gpsalt", "u", "v"]
 
             self.set_qc_ds(
@@ -246,6 +246,10 @@ def near_surface_coverage(
             )
 
         for variable in self.qc_vars.keys():
+            if variable in ["u", "v"]:
+                alt_dim = "gpsalt"
+            else:
+                alt_dim = "alt"
             dataset = ds.where(
                 (ds[alt_dim] > alt_bounds[0]) & (ds[alt_dim] < alt_bounds[1]), drop=True
             )
@@ -311,24 +315,23 @@ def low_physics(
         -------
         None
         """
-        ds = self.qc_ds
-        ds_check = ds.where(ds[alt_dim] < 100, drop=True)
+        ds_check = self.qc_ds.sortby("time", ascending=False).dropna(dim="time")
         if ds_check.sizes["time"] == 0:
             self.qc_flags["rh_low_physics"] = False
             self.qc_flags["ta_low_physics"] = False
             self.qc_flags["p_low_physics"] = False
-            self.qc_details["rh_low_physics_min"] = np.nan
-            self.qc_details["ta_low_physics_min"] = np.nan
+            self.qc_details["rh_low_physics_sfc"] = np.nan
+            self.qc_details["ta_low_physics_sfc"] = np.nan
             self.qc_details["p_low_physics_sfc"] = np.nan
         else:
-            sfc_p = ds.p.sortby("time", ascending=False).dropna(dim="time").values[0]
-            self.qc_flags["ta_low_physics"] = ds_check.ta.min() > float(ta_min)
-            self.qc_flags["rh_low_physics"] = ds_check.rh.min() > float(rh_min)
+            sfc_p = ds_check.p.values[0]
+            self.qc_flags["ta_low_physics"] = ds_check.ta.values[0] > float(ta_min)
+            self.qc_flags["rh_low_physics"] = ds_check.rh.values[0] > float(rh_min)
             self.qc_flags["p_low_physics"] = (sfc_p > float(p_min)) and (
                 sfc_p < float(p_max)
             )
-            self.qc_details["rh_low_physics_min"] = ds_check.rh.min().values
-            self.qc_details["ta_low_physics_min"] = ds_check.ta.min().values
+            self.qc_details["rh_low_physics_sfc"] = ds_check.rh.values[0]
+            self.qc_details["ta_low_physics_sfc"] = ds_check.ta.values[0]
             self.qc_details["p_low_physics_sfc"] = sfc_p
 
     def check_qc(self, used_flags=None, check_ugly=True):
@@ -547,7 +550,7 @@ def add_alt_near_gpsalt_to_ds(self, ds):
             )
 
             ds = hx.add_ancillary_var(
-                ds, "alt", "alt_near_gpsalt alt_near_gpsalt_max_diff"
+                ds, "altitude", "alt_near_gpsalt alt_near_gpsalt_max_diff"
             )
         return ds
 
@@ -571,41 +574,10 @@ def add_below_aircraft_to_ds(self, ds):
             )
         )
 
-        ds = hx.add_ancillary_var(ds, alt_dim, f"{alt_dim}_below_aircraft")
+        ds = hx.add_ancillary_var(ds, self.alt_dim, f"{alt_dim}_below_aircraft")
         return ds
 
-    def replace_alt_var(self, ds, alt_var):
-        """
-        Replace the altitude variable in a dataset with its counterpart.
-
-        This method swaps the values of the specified altitude variable with its counterpart
-        in the dataset. If `alt_var` is "alt", it will be replaced with "gpsalt", and vice versa.
-        If `alt_var` is neither "alt" nor "gpsalt", a ValueError is raised.
-
-        Parameters:
-        - ds: The dataset containing the altitude variables.
-        - alt_var: A string specifying the altitude variable to be replaced.
-                It must be either "alt" or "gpsalt".
-
-        Returns:
-        - A new dataset with the specified altitude variable replaced by its counterpart.
-
-        Raises:
-        - ValueError: If `alt_var` is not "alt" or "gpsalt".
-        """
-        if alt_var == "alt":
-            replace_var = "gpsalt"
-        elif alt_var == "gpsalt":
-            replace_var = "alt"
-        else:
-            raise ValueError(f"{alt_var} is no known altitude variable.")
-
-        ds_out = ds.assign({alt_var: ds[replace_var]})
-        self.qc_flags.update({f"{alt_var}_values": False})
-
-        return ds_out
-
-    def add_replace_alt_var_to_ds(self, ds):
+    def add_alt_source_to_ds(self, ds):
         """
         Adds  an ancillary variable in the dataset for the altitude dimension.
 
@@ -622,21 +594,16 @@ def add_replace_alt_var_to_ds(self, ds):
         - The updated dataset with the ancillary variable added or replaced.
         """
         ds = ds.assign(
-            {
-                f"{self.alt_dim}_values": np.byte(
-                    (not self.qc_flags.get(f"{self.alt_dim}_values", True))
-                )
-            }
+            {f"{self.alt_dim}_source": self.qc_flags.get(f"{self.alt_dim}_source")}
         )
-        ds[f"{self.alt_dim}_values"].attrs.update(
+        ds[f"{self.alt_dim}_source"].attrs.update(
             dict(
-                long_name=f"Values for {self.alt_dim} are present in raw data",
-                flag_values="0 1 ",
-                flag_meaning="GOOD BAD",
+                long_name=f"raw data dimension {self.alt_dim} is derived from",
+                flag_values="alt gpsalt",
             )
         )
 
-        ds = hx.add_ancillary_var(ds, self.alt_dim, f"{self.alt_dim}_values")
+        ds = hx.add_ancillary_var(ds, self.alt_dim, f"{self.alt_dim}_source")
         return ds
 
     def add_non_var_qc_to_ds(self, ds):
@@ -645,7 +612,7 @@ def add_non_var_qc_to_ds(self, ds):
 
         This method performs the following operations on the input dataset `ds`:
         1. Adds altitude near GPS altitude to the dataset using the `add_alt_near_gpsalt_to_ds` method.
-        2. Replaces altitude variable in the dataset using the `add_replace_alt_var_to_ds` method.
+        2. Replaces altitude variable in the dataset using the `add_alt_source_to_ds` method.
 
         Parameters:
         - ds: The input dataset to which non-variable QC data will be added.
@@ -654,7 +621,7 @@ def add_non_var_qc_to_ds(self, ds):
         - ds_out: The output dataset with added non-variable QC data.
         """
         ds_out = self.add_alt_near_gpsalt_to_ds(ds)
-        ds_out = self.add_replace_alt_var_to_ds(ds_out)
+        ds_out = self.add_alt_source_to_ds(ds_out)
         ds_out = self.add_below_aircraft_to_ds(ds_out)
 
         return ds_out

pydropsonde/pipeline.py

@@ -517,8 +517,8 @@ def run_pipeline(pipeline: dict, config: configparser.ConfigParser):
             "get_l2_variables",
             "convert_to_si",
             "below_aircraft_qc",
-            "replace_alt_dim",
             "get_qc",
+            "replace_alt_dim",
             #            "remove_non_qc_sondes",
         ],
         "output": "sondes",

pydropsonde/processor.py

@@ -1041,7 +1041,7 @@ def remove_above_aircraft(self, max_alt=15000):
         """
         remove measured values above aircraft
         """
-        variables = ["lat", "lon", "gpsalt", "u", "v"]
+        variables = ["lat", "lon", self.alt_dim, "u", "v"]
         maxalt = self.flight_attrs.get("aircraft_msl_altitude_(m)", float(max_alt))
         self.interim_l3_ds = hx.remove_above_alt(
             self.interim_l3_ds, variables, alt_dim=self.alt_dim, maxalt=maxalt
@@ -1103,7 +1103,9 @@ def add_iwv(self):
         self : object
             Returns the sonde object with integrated water vapour added to the interim l3 dataset.
         """
-        self.interim_l3_ds = hh.calc_iwv(self.interim_l3_ds, qc_var=["rh_qc", "ta_qc"])
+        self.interim_l3_ds = hh.calc_iwv(
+            self.interim_l3_ds, qc_var=["rh_qc", "ta_qc"], alt_dim=self.alt_dim
+        )
 
         return self
 
@@ -1160,42 +1162,52 @@ def set_alt_dim(self, alt_dim="alt"):
 
     def replace_alt_dim(self, drop_nan=True):
         """
-        Replaces the altitude dimension in the dataset if all altitude values are NaN.
+        Replaces the altitude dimension in the dataset if one altitude coordinate is worse than the other
 
-        This function loads the interim level 2 dataset and checks the altitude variable for NaN values.
-        If all values are NaN, it attempts to replace the altitude variable with the other (alt and gpsalt)
-        Depending on the `drop_nan` parameter, it either drops sonde from the dictionary or sets the qc values.
-        If the altitude values remain NaN after replacement, the dataset is dropped.
+        -> if no gpsalt values are present and the sonde probably didn't reach the ground it is dropped
+        -> if the sonde did not reach the ground, but gpsalt is available, gpsalt is used
+        -> if gpsalt does not go until the surface, but pressure measurements suggest that the sonde reached the surface, alt is used
 
-        Parameters:
-        - drop_nan (bool): Determines whether to drop the dataset if altitude values are NaN.
-        If True, the dataset is dropped; otherwise, it attempts to switch the altitude variable.
 
         Returns:
         - self: Returns the object itself if the altitude dimension is successfully replaced or remains valid.
         - None: Returns None if the dataset is dropped due to NaN altitude values.
         """
         alt_dim = self.alt_dim
-        ds = self.interim_l2_ds.load()
-        alt = ds[alt_dim]
-        if np.all(np.isnan(alt)):
-            ds = self.qc.replace_alt_var(ds, alt_dim)
-            if hh.get_bool(drop_nan):
-                print(
-                    f"No {alt_dim} values.  Sonde {self.serial_id} from {self.flight_id} is dropped"
-                )
-                return None
-            else:
-                print(
-                    f"No {alt_dim} values.  Sonde {self.serial_id} from {self.flight_id} altitude is switched"
-                )
-                alt = ds[alt_dim]
-                if np.all(np.isnan(alt)):
-                    print(
-                        f"No altitude values. Sonde {self.serial_id} from {self.flight_id} is dropped"
-                    )
-                    return None
-                self.interim_l2_ds = ds
+        ds = self.interim_l2_ds
+        if (not self.qc.qc_flags["p_low_physics"]) and (np.all(np.isnan(ds["gpsalt"]))):
+            print(
+                f"No gpsalt values and no reliable alt values.  Sonde {self.serial_id} from {self.flight_id} is dropped"
+            )
+            return None
+        elif alt_dim == "alt":
+            self.qc.qc_flags.update({"altitude_source": "alt"})
+            if not self.qc.qc_flags["p_low_physics"]:
+                for var in ["rh", "ta", "p"]:
+                    self.qc.qc_flags[f"{var}_near_surface"] = False
+                    self.qc.qc_details[f"{var}_near_surface_count"] = np.nan
+
+                ds = ds.assign({"alt": ds["gpsalt"]})
+                self.qc.qc_flags.update({"altitude_source": "gpsalt"})
+            ds = ds.rename({"alt": "altitude"}).drop_vars(["gpsalt"])
+
+        elif alt_dim == "gpsalt":
+            self.qc.qc_flags.update({"altitude_source": "gpsalt"})
+            if (not self.qc.qc_flags["u_near_surface"]) and (
+                self.qc.qc_flags["p_low_physics"]
+            ):
+                ds = ds.assign({alt_dim: ds["alt"]})
+                self.qc.qc_flags.update({"altitude_source": "alt"})
+            elif not self.qc.qc_flags["p_low_physics"]:
+                for var in ["rh", "ta", "p"]:
+                    self.qc.qc_flags[f"{var}_near_surface"] = False
+                    self.qc.qc_details[f"{var}_near_surface_count"] = np.nan
+            ds = ds.rename({"gpsalt": "altitude"}).drop_vars(["alt"])
+        else:
+            self.qc.qc_flags.update({"altitude_source": self.alt_dim})
+        self.alt_dim = "altitude"
+        self.qc.alt_dim = "altitude"
+        self.interim_l2_ds = ds
         return self
 
     def swap_alt_dimension(self):
@@ -1609,7 +1621,7 @@ def add_qc_to_interim_l3(self, keep=None):
                 keep = (
                     [f"{var}_qc" for var in list(self.qc.qc_by_var.keys())]
                     + list(self.qc.qc_details.keys())
-                    + ["alt_near_gpsalt"]
+                    + ["alt_near_gpsalt", "altitude_source"]
                 )
                 for variable in self.qc.qc_vars:
                     ds = self.qc.add_variable_flags_to_ds(ds, variable, details=True)
@@ -1625,6 +1637,7 @@ def add_qc_to_interim_l3(self, keep=None):
                     ds = self.qc.add_variable_flags_to_ds(
                         ds, "ta", add_to="theta", details=True
                     )
+                ds = self.qc.add_non_var_qc_to_ds(ds)
             elif keep == "var_flags":
                 keep = [f"{var}_qc" for var in list(self.qc.qc_by_var.keys())] + [
                     "sonde_qc"