Feature/individual epsilons (#29)

* Add individual epsilons for obs types for OI

* Feature/cryoclim (#28)

* Work by Jostein on cryoclim in pysurfex

* No fixed version

* Add cycolim/sentinel as an observation set

* Correct linting

* remove prints

* Remove comma

* Add option to read input provided variable name in cryoclim file

* Reformat

* Install in local environment with poetry (#30)

* Avoid local test files

* Bad merge

* Enable sigmao also for clients. Add sigmao also to json obs sets.

* Sigmaos must also enter the observation objects in the QC dataset

pysurfex/bufr.py

@@ -36,6 +36,7 @@ def __init__(
         latrange=None,
         label="bufr",
         use_first=False,
+        sigmao=None,
     ):
         """Initialize a bufr observation set.
 
@@ -48,6 +49,7 @@ def __init__(
             latrange (list): Allowed range of latitides [min, max]
             label (str): A label for the resulting observations set
             use_first (bool): Use only the first valid observation for a point if more are found
+            sigmao (float, optional): Observation error relative to normal background error. Defaults to None.
 
         Raises:
             RuntimeError: ECCODES not found. Needed for bufr reading
@@ -438,7 +440,7 @@ def __init__(
         # close the file
         file_handler.close()
 
-        ObservationSet.__init__(self, observations, label=label)
+        ObservationSet.__init__(self, observations, label=label, sigmao=sigmao)
 
     @staticmethod
     def td2rh(t_d, temp, kelvin=True):

pysurfex/cli.py

@@ -799,7 +799,7 @@ def run_gridpp(**kwargs):
         elev_gradient = kwargs["elev_gradient"]
     if elev_gradient != -0.0065 and elev_gradient != 0:
         raise RuntimeError("Not a valid elevation gradient")
-    epsilon = 0.25
+    epsilon = None
     if "epsilon" in kwargs:
         epsilon = kwargs["epsilon"]
     minvalue = None
@@ -822,7 +822,6 @@ def run_gridpp(**kwargs):
     # Read OK observations
     observations = dataset_from_file(an_time, obs_file, qc_flag=0)
     logging.info("Found %s observations with QC flag == 0", str(len(observations.lons)))
-
     field = horizontal_oi(
         geo,
         background,
@@ -1812,6 +1811,7 @@ def bufr2json(argv=None):
     output = kwargs.get("output")
     valid_dtg = as_datetime(kwargs.get("dtg"))
     valid_range = as_timedelta(seconds=kwargs.get("valid_range"))
+    sigmao = kwargs.get("sigmao")
     label = kwargs.get("label")
     indent = kwargs.get("indent")
     lonrange = kwargs.get("lonrange")
@@ -1828,6 +1828,7 @@ def bufr2json(argv=None):
         latrange=latrange,
         label=label,
         indent=indent,
+        sigmao=sigmao,
     )
 
 
@@ -1862,6 +1863,7 @@ def obs2json(argv=None):
     level = kwargs.get("level")
     obtypes = kwargs.get("obtypes")
     subtypes = kwargs.get("subtypes")
+    sigmao = kwargs.get("sigmao")
     pos_t_range = kwargs.get("pos_t_range")
     neg_t_range = kwargs.get("neg_t_range")
     indent = kwargs.get("indent")
@@ -1888,6 +1890,7 @@ def obs2json(argv=None):
         indent=indent,
         obtypes=obtypes,
         subtypes=subtypes,
+        sigmao=sigmao,
     )
 
 

pysurfex/cmd_parsing.py

@@ -1079,7 +1079,7 @@ def parse_args_gridpp(argv):
         choices=[0, -0.0065],
     )
     parser.add_argument(
-        "--epsilon", dest="epsilon", type=float, default=0.25, required=False
+        "--epsilon", dest="epsilon", type=float, default=None, required=False
     )
     parser.add_argument(
         "--minvalue", dest="minvalue", type=float, default=None, required=False
@@ -1387,6 +1387,14 @@ def parse_args_bufr2json(argv):
         help="Valid range in seconds",
         default=3600,
     )
+    parser.add_argument(
+        "--sigmao",
+        dest="sigmao",
+        type=float,
+        required=False,
+        default=None,
+        help="Observation error relative to normal background error.",
+    )
     parser.add_argument(
         "--debug", action="store_true", help="Debug", required=False, default=False
     )
@@ -1484,6 +1492,14 @@ def parse_args_obs2json(argv):
         default=None,
         help="Subtypes (obsoul)",
     )
+    parser.add_argument(
+        "--sigmao",
+        dest="sigmao",
+        type=float,
+        required=False,
+        default=None,
+        help="Observation error relative to normal background error.",
+    )
     parser.add_argument(
         "--debug", action="store_true", help="Debug", required=False, default=False
     )
@@ -2033,7 +2049,7 @@ def parse_cryoclim_pseudoobs(argv):
         required=False,
     )
     parser.add_argument(
-        "-gt",
+        "-lt",
         "--laf_threshold",
         dest="laf_threshold",
         type=float,

pysurfex/input_methods.py

@@ -59,6 +59,8 @@ def get_datasources(obs_time, settings):
                     kwargs.update({"lonrange": settings[obs_set]["lonrange"]})
                 if "latrange" in settings[obs_set]:
                     kwargs.update({"latrange": settings[obs_set]["latrange"]})
+                if "sigmao" in settings[obs_set]:
+                    kwargs.update({"sigmao": settings[obs_set]["sigmao"]})
                 if "dt" in settings[obs_set]:
                     deltat = settings[obs_set]["dt"]
                 else:
@@ -87,6 +89,8 @@ def get_datasources(obs_time, settings):
                         pos_t_range = 15
                         if "pos_t_range" in settings[obs_set]:
                             pos_t_range = settings[obs_set]["pos_t_range"]
+                        if "sigmao" in settings[obs_set]:
+                            kwargs.update({"sigmao": settings[obs_set]["sigmao"]})
 
                         dtg = validtime - as_timedelta(seconds=int(neg_t_range) * 60)
                         end_dtg = validtime + as_timedelta(seconds=int(pos_t_range) * 60)
@@ -111,6 +115,8 @@ def get_datasources(obs_time, settings):
                     kwargs.update({"lonrange": settings[obs_set]["lonrange"]})
                 if "latrange" in settings[obs_set]:
                     kwargs.update({"latrange": settings[obs_set]["latrange"]})
+                if "sigmao" in settings[obs_set]:
+                    kwargs.update({"sigmao": settings[obs_set]["sigmao"]})
                 if "dt" in settings[obs_set]:
                     kwargs.update({"dt": settings[obs_set]["dt"]})
                 else:
@@ -135,6 +141,8 @@ def get_datasources(obs_time, settings):
                     kwargs.update({"latrange": settings[obs_set]["latrange"]})
                 if "unit" in settings[obs_set]:
                     kwargs.update({"unit": settings[obs_set]["unit"]})
+                if "sigmao" in settings[obs_set]:
+                    kwargs.update({"sigmao": settings[obs_set]["sigmao"]})
                 if "level" in settings[obs_set]:
                     kwargs.update({"level": settings[obs_set]["level"]})
                 kwargs.update({"validtime": obs_time})
@@ -150,6 +158,7 @@ def get_datasources(obs_time, settings):
                 pos_dt = None
                 obtypes = None
                 subtypes = None
+                sigmao = None
                 if "obnumber" in settings[obs_set]:
                     obnumber = int(settings[obs_set]["obnumber"])
                 if "neg_dt" in settings[obs_set]:
@@ -160,6 +169,8 @@ def get_datasources(obs_time, settings):
                     obtypes = settings[obs_set]["obtypes"]
                 if "subtypes" in settings[obs_set]:
                     subtypes = settings[obs_set]["subtypes"]
+                if "sigmao" in settings[obs_set]:
+                    sigmao = settings[obs_set]["sigmao"]
                 if os.path.exists(filename):
                     datasources.append(
                         ObservationDataSetFromObsoulFile(
@@ -170,6 +181,7 @@ def get_datasources(obs_time, settings):
                             obtypes=obtypes,
                             subtypes=subtypes,
                             obnumber=obnumber,
+                            sigmao=sigmao,
                         )
                     )
                 else:
@@ -185,6 +197,8 @@ def get_datasources(obs_time, settings):
                     varname = settings[obs_set]["varname"]
 
                 kwargs.update({"var": varname})
+                if "sigmao" in settings[obs_set]:
+                    kwargs.update({"sigmao": settings[obs_set]["sigmao"]})
                 if os.path.exists(filename):
                     datasources.append(JsonObservationSet(filename, **kwargs))
                 else:
@@ -230,7 +244,7 @@ def set_geo_from_obs_set(
 
     logging.debug("%s", settings)
     logging.debug("Get data source")
-    __, lons, lats, __, __, __, __ = get_datasources(obs_time, settings)[0].get_obs()
+    __, lons, lats, __, __, __, __, __ = get_datasources(obs_time, settings)[0].get_obs()
 
     selected_lons = []
     selected_lats = []
@@ -271,6 +285,7 @@ def get_obsset(
     level=None,
     obtypes=None,
     subtypes=None,
+    sigmao=None,
 ):
     """Create an observation set from an input data set.
 
@@ -288,6 +303,7 @@ def get_obsset(
         level (str, optional): Level (FROST)
         obtypes (list, optional): Obstypes (obsoul)
         subtypes (list, optional): Subtypes (obsoul)
+        sigmao (float, optional): Observation error relative to normal background error. Defaults to None.
 
     Returns:
         obsset (ObservationSet): Observation set
@@ -331,6 +347,7 @@ def get_obsset(
             "level": level,
             "obtypes": obtypes,
             "subtypes": subtypes,
+            "sigmao": sigmao,
             "pos_t_range": pos_t_range_seconds,
             "neg_t_range": neg_t_range_seconds,
         }
@@ -357,6 +374,7 @@ def create_obsset_file(
     level=None,
     obtypes=None,
     subtypes=None,
+    sigmao=None,
 ):
     """Create an observation set from an input data set.
 
@@ -371,11 +389,12 @@ def create_obsset_file(
         lonrange (tuple, optional): Longitude range (min, max). Defaults to None.
         latrange (tuple, optional): Latitude range (min, max). Defaults to None.
         label (str, optional): Obs set label. Default to None which means it will be the same as obs_type
-        indent (int, optional): File indentation. Default to None.
+        indent (int, optional): File indentation. Defaults to None.
         unit (str, optional): Unit (FROST)
         level (str, optional): Level (FROST)
         obtypes (list, optional): Obstypes (obsoul)
         subtypes (list, optional): Subtypes (obsoul)
+        sigmao (float, optional): Observation error relative to normal background error. Defaults to None.
 
     """
     logging.debug("Get data source")
@@ -393,5 +412,6 @@ def create_obsset_file(
         level=level,
         obtypes=obtypes,
         subtypes=subtypes,
+        sigmao=sigmao,
     )
     obsset.write_json_file(output, indent=indent)

pysurfex/interpolation.py

@@ -455,7 +455,7 @@ def horizontal_oi(
     elev_gradient=None,
     structure_function="Barnes",
     max_locations=50,
-    epsilon=0.5,
+    epsilon=None,
     minvalue=None,
     maxvalue=None,
     interpol="bilinear",
@@ -503,12 +503,12 @@ def obs2vectors(my_obs):
             my_obs.stids,
             my_obs.elevs,
             my_obs.values,
-            my_obs.cis,
+            my_obs.epsilons,
             my_obs.lafs,
         )
 
     vectors = np.vectorize(obs2vectors)
-    lons, lats, __, elevs, values, __, __ = vectors(observations)
+    lons, lats, __, elevs, values, sigmaos, __ = vectors(observations)
 
     bgrid = Grid(glons, glats, gelevs)
     points = Points(lons, lats, elevs)
@@ -523,26 +523,36 @@ def obs2vectors(my_obs):
         lats2 = []
         elevs2 = []
         values2 = []
+        sigmaos2 = []
         for point in range(0, len(lons)):
             if np.isnan(pbackground[point]):
                 logging.info(
                     "Undefined background in lon=%s lat=%s value=%s",
                     lons[point],
                     lats[point],
                     values[point],
+                    sigmaos[point],
                 )
             else:
                 lons2.append(lons[point])
                 lats2.append(lats[point])
                 elevs2.append(elevs[point])
                 values2.append(values[point])
+                sigmaos2.append(sigmaos[point])
         values = values2
+        sigmaos = sigmaos2
         points = Points(lons2, lats2, elevs2)
         pbackground = grid2points(
             bgrid, points, background, operator=interpol, elev_gradient=elev_gradient
         )
 
-    variance_ratios = np.full(points.points.size(), epsilon)
+    # Set relationship between obs/background error
+    if epsilon is None:
+        logging.info("Using epsilon from observation data sets")
+        variance_ratios = sigmaos
+    else:
+        logging.info("Using fixed epsilon %s", epsilon)
+        variance_ratios = np.full(points.points.size(), epsilon)
 
     if structure_function == "Barnes":
         structure = gridpp.BarnesStructure(hlength, vlength, wlength)