fixed flags issues for water after pull master branch

Author: clemgoyens

hypernets_processor/calibration/calibrate.py

@@ -74,6 +74,7 @@ def calibrate_l1a(self, measurandstring, dataset_l0, dataset_l0_bla):
     def prepare_calibration_data(self,measurandstring):
         hypstar=self.context.get_config_value("hypstar_cal_number")
         directory=self.context.get_config_value("calibration_directory")
+        print(directory)
         caldates=[os.path.basename(path) for path in glob.glob
         (os.path.join(directory,"hypstar_"+str(hypstar)+"/radiometric/*"))]
         caldate=caldates[0]

hypernets_processor/etc/job_CG.config

@@ -0,0 +1,52 @@
+[Job]
+job_name = run_20201028T150017
+job_working_directory = home/cgoyens/.hypernets/tmp
+
+#SEQ20201016T141413
+sequence_id: SEQ20200821T154756
+Lat: 43.69886
+Lon: 7.30720
+site_abbr: test
+
+[Input]
+raw_data_directory = /home/cgoyens/OneDrive/BackUpThinkpadClem/Projects/HYPERNETS/NetworkDesign_D52/DataProcChain/hypernets_processor/hypernets_processor/data_io/tests/reader/SEQ20200821T154756
+[Log]
+log_path = test.log
+verbose = True
+quiet = False
+
+[Processor]
+network = l
+
+[Output]
+product_format = netcdf
+write_l0 = False
+write_l1a = False
+write_l1b = True
+write_l1c = False
+write_l2a = True
+archive_directory = home/cgoyens/OneDrive/BackUpThinkpadClem/Projects/HYPERNETS/NetworkDesign_D52/DataProcChain/hypernets_processor/hypernets_processor/test_hypernets
+
+
+[Plotting]
+plotting_directory: ./out/plots/
+plotting_format: png
+plot_l1a: True
+plot_l1b: True
+plot_l2a: True
+
+[Calibration]
+calibration_directory: /home/cgoyens/OneDrive/BackUpThinkpadClem/Projects/HYPERNETS/NetworkDesign_D52/DataProcChain/hypernets_processor/calibration_files/HYPSTAR_cal
+hypstar_cal_number:220241
+cal_date:200728
+measurement_function_calibrate: StandardMeasurementFunction
+
+[Combine VIS-SWIR]
+measurement_function_combine: StepCombine
+wavelength_switch: 1044
+
+[Interpolate]
+measurement_function_interpolate: WaterNetworkInterpolationLinear
+
+[SurfaceReflectance]
+measurement_function_surface_reflectance = WaterNetworkProtocol

hypernets_processor/rhymer/rhymer/hypstar/rhymer_hypstar.py

@@ -16,6 +16,7 @@
 from hypernets_processor.plotting.plotting import Plotting
 from hypernets_processor.data_io.format.metadata import METADATA_DEFS
 from hypernets_processor.data_io.format.variables import VARIABLES_DICT_DEFS
+from hypernets_processor.data_io.dataset_util import DatasetUtil as du
 
 import numpy as np
 import math
@@ -29,7 +30,7 @@ def __init__(self, context):
         self.writer = HypernetsWriter(context)
         self.cal = Calibrate(context, MCsteps=100)
         self.intp = Interpolate(context, MCsteps=1000)
-        self.plot=Plotting(context)
+        self.plot = Plotting(context)
         self.rhymeranc = RhymerAncillary(context)
         self.rhymerproc = RhymerProcessing(context)
         self.rhymershared = RhymerShared(context)
@@ -89,26 +90,18 @@ def qc_scan(self, dataset, measurandstring):
                 if v > self.context.get_config_value("diff_threshold"):
                     # get flag value for the temporal variability
                     if measurandstring == 'irradiance':
-                        flagval=2 ** (self.context.get_config_value("temp_variability_ed"))
+                        dataset['quality_flag'][scans[i]] = du.set_flag(dataset["quality_flag"][scans[i]],
+                                                                        "temp_variability_ed")
                     else:
-                        flagval = 2 ** (self.context.get_config_value("temp_variability_lu"))
+                        dataset['quality_flag'][scans[i]] = du.set_flag(dataset["quality_flag"][scans[i]],
+                                                                        "temp_variability_lu")
 
-                    dataset['quality_flag'][scans[i]] = dataset['quality_flag'][scans[i]]+flagval
                     seq = dataset.attrs["sequence_id"]
                     ts = datetime.utcfromtimestamp(dataset['acquisition_time'][i])
 
-                    if verbosity > 2: self.context.logger.info('Temporal jump: in {}:  Aquisition time {}, {}'.format(seq, ts, ', '.join(
-                        ['{}:{}'.format(k, dataset[k][scans[i]].values) for k in ['scan', 'quality_flag']])))
-
-                # ## check for complete spectra
-                # ## assumed here they are masked by inf or nan by calibration step
-                # if any([not np.isfinite(j) for j in data[i]]):
-                #     # set quality flag to 10
-                #     flagval = 2 ** (self.context.get_config_value("inf_nan"))
-                #     dataset['quality_flag'][scans[i]] = dataset['quality_flag'][scans[i]] + flagval
-                #     ts = datetime.utcfromtimestamp(dataset['acquisition_time'][i])
-                #     if verbosity > 2: print('inf or Nan data in: Aquisition time {}, {}'.format(ts, ', '.join(
-                #         ['{}:{}'.format(k, dataset[k][scans[i]].values) for k in ['scan', 'quality_flag']])))
+                    if verbosity > 2: self.context.logger.info(
+                        'Temporal jump: in {}:  Aquisition time {}, {}'.format(seq, ts, ', '.join(
+                            ['{}:{}'.format(k, dataset[k][scans[i]].values) for k in ['scan', 'quality_flag']])))
 
             return dataset
 
@@ -119,9 +112,6 @@ def cycleparse(self, rad, irr):
         nbrlu = self.context.get_config_value("n_upwelling_rad")
         nbred = self.context.get_config_value("n_upwelling_irr")
         nbrlsky = self.context.get_config_value("n_downwelling_rad")
-        flag_nbrlu = 2**self.context.get_config_value("min_nbrlu")
-        flag_nbred = 2**self.context.get_config_value("min_nbred")
-        flag_nbrlsky = 2**self.context.get_config_value("min_nbrlsky")
 
         if protocol != 'WaterNetworkProtocol':
             # here we should simply provide surface reflectance?
@@ -142,12 +132,9 @@ def cycleparse(self, rad, irr):
                     senz = float(senz)
                     sena = abs(float(sena))
                 else:
-                    flagval = 2 ** (self.context.get_config_value("angles_missing"))
-                    scani['quality_flag'] = scani['quality_flag'] + flagval
-                    if self.context.get_config_value("verbosity") > 2:
-                        self.context.logger.info(
-                        'NULL angles: Aquisition time {}, {}'.format(ts, ', '.join(
-                            ['{}:{}'.format(k, scani[k].values) for k in ['scan', 'quality_flag']])))
+                    scani['quality_flag'] = du.set_flag(scani['quality_flag'], "angles_missing")
+                    self.context.logger.info('NULL angles: Aquisition time {}, {}'.format(ts, ', '.join(
+                        ['{}:{}'.format(k, scani[k].values) for k in ['scan', 'quality_flag']])))
                     continue
 
                 # ## identify measurement
@@ -169,46 +156,47 @@ def cycleparse(self, rad, irr):
             sena_lsky = np.unique(lsky["viewing_azimuth_angle"].values)
             for i in sena_lu:
                 if i not in sena_lsky:
-                    flagval=self.context.get_config_value("lu_eq_missing")
-                    lu["quality_flag"][lu["viewing_azimuth_angle"] == i] = lu["quality_flag"][lu["viewing_azimuth_angle"] == i]+2**flagval
+                    lu["quality_flag"][lu["viewing_azimuth_angle"] == i] = du.set_flag(
+                        lu["quality_flag"][lu["viewing_azimuth_angle"] == i], "lu_eq_missing")
                     if self.context.get_config_value("verbosity") > 2:
                         ts = [datetime.utcfromtimestamp(x) for x in
                               lu['acquisition_time'][lu["viewing_azimuth_angle"] == i].values]
-                        self.context.logger.info('No azimuthal equivalent downwelling radiance measurement: Aquisition time {}, {}'.format(
-                            ts, ', '.join(
-                                ['{}:{}'.format(k, lu[k][lu["viewing_azimuth_angle"] == i].values) for k in
-                                 ['scan', 'quality_flag']])))
+                        self.context.logger.info(
+                            'No azimuthal equivalent downwelling radiance measurement: Aquisition time {}, {}'.format(
+                                ts, ', '.join(
+                                    ['{}:{}'.format(k, lu[k][lu["viewing_azimuth_angle"] == i].values) for k in
+                                     ['scan', 'quality_flag']])))
 
             # check if we have the required fresnel angle for lsky
             senz_lu = np.unique(lu["viewing_zenith_angle"].values)
             senz_lsky = 180 - np.unique(lsky["viewing_zenith_angle"].values)
             for i in senz_lu:
                 if i not in senz_lsky:
-                    flagval = self.context.get_config_value("fresnel_angle_missing")
-                    lu["quality_flag"][lu["viewing_azimuth_angle"] == i] = lu["quality_flag"][lu["viewing_azimuth_angle"] == i] + 2 ** flagval
-                    if self.context.get_config_value("verbosity") > 2:
-                        ts = [datetime.utcfromtimestamp(x) for x in
-                              lu['acquisition_time'][lu["viewing_zenith_angle"] == i].values]
-                        self.context.logger.info(
+                    lu["quality_flag"][lu["viewing_azimuth_angle"] == i] = du.set_flag(lu["quality_flag"][lu["viewing_azimuth_angle"]==i],"fresnel_angle_missing")
+                    ts = [datetime.utcfromtimestamp(x) for x in lu['acquisition_time'][lu["viewing_zenith_angle"] == i].values]
+                    self.context.logger.info(
                             'No downwelling radiance measurement at appropriate fresnel angle: Aquisition time {}, {}'.format(
                                 ts, ', '.join(
                                     ['{}:{}'.format(k, lu[k][lu["viewing_azimuth_angle"] == i].values) for k
                                      in ['scan', 'quality_flag']])))
 
             # check if correct number of radiance and irradiance data
+
             if lu.scan[lu['quality_flag'] <= 0].count() < nbrlu:
-                lu["quality_flag"].values = [lu.sel(scan=i)["quality_flag"].values+flag_nbrlu for i in lu['scan']]
-                if self.context.get_config_value("verbosity") > 2:
-                    self.context.logger.info("No enough upwelling radiance data for sequence {}".format(lu.attrs['sequence_id']))
+                lu["quality_flag"].values = [du.set_flag(lu["quality_flag"][lu["scan"] == i], "min_nbrlu") for
+                                             i in lu["scan"]]
+                self.context.logger.info(
+                    "No enough upwelling radiance data for sequence {}".format(lu.attrs['sequence_id']))
             if lsky.scan[lsky['quality_flag'] <= 1].count() < nbrlsky:
-                lsky["quality_flag"].values = [lsky.sel(scan=i)["quality_flag"].values+flag_nbrlsky for i in lsky['scan']]
-                if self.context.get_config_value("verbosity") > 2:
-                    self.context.logger.info("No enough downwelling radiance data for sequence {}".format(lsky.attrs['sequence_id']))
-
+                lsky["quality_flag"].values = [du.set_flag(lsky["quality_flag"][lu["scan"] == i], "min_nbrlu") for
+                                               i in lsky["scan"]]
+                self.context.logger.info(
+                    "No enough downwelling radiance data for sequence {}".format(lsky.attrs['sequence_id']))
             if irr.scan[irr['quality_flag'] <= 1].count() < nbred:
-                irr["quality_flag"].values = [irr.sel(scan=i)["quality_flag"].values + 2 ** flag_nbred for i in irr['scan']]
-                if self.context.get_config_value("verbosity") > 2:
-                    self.context.logger.info("No enough irradiance data for sequence {}".format(irr.attrs['sequence_id']))
+                irr["quality_flag"].values = [du.set_flag(lsky["quality_flag"][irr["scan"] == i], "min_nbrlu") for
+                                              i in irr["scan"]]
+                self.context.logger.info(
+                    "No enough downwelling irradiance data for sequence {}".format(irr.attrs['sequence_id']))
 
             return lu, lsky, irr
 
@@ -217,11 +205,11 @@ def get_wind(self, l1b):
         lat = l1b.attrs['site_latitude']
         lon = l1b.attrs['site_latitude']
         wind = []
-        flagval = self.context.get_config_value("def_wind_flag")
         for i in range(len(l1b.scan)):
             wa = self.context.get_config_value("wind_ancillary")
-            if wa == False:
-                l1b.sel(scan=i)["quality_flag"]=l1b.sel(scan=i)["quality_flag"] + 2 ** flagval
+            if not wa:
+                l1b["quality_flag"][l1b["scan"] == i] = du.set_flag(l1b["quality_flag"][l1b["scan"] == i],
+                                                                    "def_wind_flag")
                 self.context.logger.info("Default wind speed {}".format(self.context.get_config_value("wind_default")))
                 wind.append(self.context.get_config_value("wind_default"))
             else:
@@ -261,18 +249,16 @@ def get_fresnelrefl(self, l1b):
                 if (fresnel_sza[i] is not None) & (fresnel_raa[i] is not None):
                     sza = min(fresnel_sza[i], 79.999)
                     rhof = self.rhymerproc.mobley_lut_interp(sza, fresnel_vza[i], fresnel_raa[i],
-                                                                wind=wind[i])
+                                                             wind=wind[i])
                 else:
-                    # add a quality flag!
-                    flagval = self.context.get_config_value("fresnel_default")
-                    l1b["quality_flag"][i].values = l1b["quality_flag"][i].values + 2 ** flagval
+                    l1b["quality_flag"][l1b["scan"] == i] = du.set_flag(l1b["quality_flag"][l1b["scan"] == i],
+                                                                        "fresnel_default")
                     rhof = self.context.get_config_value("rhof_default")
-
             if self.context.get_config_value("fresnel_option") == 'Ruddick2006':
                 rhof = self.context.get_config_value("rhof_default")
                 self.context.logger.info("Apply Ruddick et al., 2006")
                 if wind[i] is not None:
-                    rhof = rhof + 0.00039 * wind[i]  + 0.000034 * wind[i]  ** 2
+                    rhof = rhof + 0.00039 * wind[i] + 0.000034 * wind[i] ** 2
 
             fresnel_coeff[i] = rhof
 
@@ -286,7 +272,7 @@ def get_fresnelrefl(self, l1b):
     def get_epsilon(self, rhow_nosc, wavelength):
 
         # wavelength = l1b['wavelength'].values
-        #get length of transposed rhow_nosc (1 epsilon per scan!)
+        # get length of transposed rhow_nosc (1 epsilon per scan!)
         epsilon = np.zeros(len(rhow_nosc.T))
         failSimil = np.zeros(len(rhow_nosc.T))
 
@@ -309,10 +295,9 @@ def get_epsilon(self, rhow_nosc, wavelength):
                 else:
                     if self.context.get_config_value("verbosity") > 2: self.context.logger.info('Passed simil test.')
             epsilon[i] = eps
-            failSimil[i]=fail_simil
+            failSimil[i] = fail_simil
         return epsilon, failSimil
 
-
     def get_rhow_nosc(self, l1b):
 
         ## read mobley rho lut
@@ -321,8 +306,8 @@ def get_rhow_nosc(self, l1b):
         fresnel_coeff = l1b['rhof'].values
         rhow_nosc_all = np.zeros((len(l1b.scan), len(wavelength)))
         lw_all = np.zeros((len(l1b.scan), len(wavelength)))
-        #rhow_all = np.zeros((len(l1b.scan), len(wavelength)))
-        #epsilon = np.zeros(len(l1b.scan))
+        # rhow_all = np.zeros((len(l1b.scan), len(wavelength)))
+        # epsilon = np.zeros(len(l1b.scan))
         simil_flag = np.zeros(len(l1b.scan))
 
         for i in range(len(l1b.scan)):
@@ -381,7 +366,7 @@ def process_l1b(self, L1a_rad, L1a_irr):
         L1b_irr = self.cal.average_l1b("irradiance", L1a_irr)
         # INTERPOLATE Lsky and Ed FOR EACH Lu SCAN! Threshold in time -> ASSIGN FLAG
         L1b = self.intp.interpolate_l1b_w(L1a_uprad, L1b_downrad, L1b_irr)
-        if self.context.get_config_value("write_l1b")==True:
+        if self.context.get_config_value("write_l1b") == True:
             self.writer.write(L1b, overwrite=True)
         return L1b
 
@@ -393,21 +378,21 @@ def process_l1c(self, l1b):
 
         rhow_nosc_all, lw_all = self.get_rhow_nosc(dataset_l1c)
 
-        #dataset_l1c['reflectance'].values = rhow_all.T
+        # dataset_l1c['reflectance'].values = rhow_all.T
 
         dataset_l1c['reflectance_nosc'].values = rhow_nosc_all.T
-        #dataset_l1c['epsilon'].values = epsilon
+        # dataset_l1c['epsilon'].values = epsilon
         dataset_l1c['water_leaving_radiance'].values = lw_all.T
 
         if self.context.get_config_value("write_l1c"):
-            self.writer.write(dataset_l1c,overwrite=True)
+            self.writer.write(dataset_l1c, overwrite=True)
 
         if self.context.get_config_value("plot_l1c"):
-            self.plot.plot_scans_in_series("irradiance",dataset_l1c)
+            self.plot.plot_scans_in_series("irradiance", dataset_l1c)
 
         return dataset_l1c
 
-    def l1c_from_l1b_dataset(self,dataset_l1b):
+    def l1c_from_l1b_dataset(self, dataset_l1b):
         """
         Makes a L2 template of the data, and propagates the appropriate keywords from L1.
 
@@ -416,10 +401,10 @@ def l1c_from_l1b_dataset(self,dataset_l1b):
         :return:
         :rtype:
         """
-        l1c_dim_sizes_dict = {"wavelength":len(dataset_l1b["wavelength"]),
-                              "scan":len(np.unique(dataset_l1b['scan']))}
+        l1c_dim_sizes_dict = {"wavelength": len(dataset_l1b["wavelength"]),
+                              "scan": len(np.unique(dataset_l1b['scan']))}
 
-        dataset_l1c = self.hdsb.create_ds_template(l1c_dim_sizes_dict,"W_L1C",
+        dataset_l1c = self.hdsb.create_ds_template(l1c_dim_sizes_dict, "W_L1C",
                                                    propagate_ds=dataset_l1b)
         dataset_l1c = dataset_l1c.assign_coords(wavelength=dataset_l1b.wavelength)
 

hypernets_processor/surface_reflectance/surface_reflectance.py

@@ -47,11 +47,9 @@ def process_l1d(self, dataset_l1c):
                                                   dataset_l1d["wavelength"].values)
             dataset_l1d["epsilon"].values = epsilon
 
-            flagval = 2 ** (self.context.get_config_value("simil_fail"))
-
-            dataset_l1d["quality_flag"].values = [
-                flagval + dataset_l1d["quality_flag"].values[i] if failSimil[i] == True else
-                dataset_l1d["quality_flag"].values[i] for i in range(len(failSimil))]
+            dataset_l1d["quality_flag"][dataset_l1d["scan"] == [i for i, x in enumerate(failSimil) if x]] = du.set_flag(
+                dataset_l1d["quality_flag"][dataset_l1d["scan"] == [i for i, x in enumerate(failSimil) if x]],
+                "simil_fail")
 
             input_vars = l1ctol1d_function.get_argument_names()
             input_qty = self.find_input(input_vars, dataset_l1d)
@@ -260,8 +258,13 @@ def l2_from_l1d_dataset(self, datasetl1d):
                                    "solar_zenith_angle"]:
                 temp_arr = np.empty(len(series_id))
                 for i in range(len(series_id)):
-                    ids = np.where((datasetl1d['series_id'] == series_id[i]) & (
-                            datasetl1d['quality_flag'] == 1))
+                    flags = ["saturation", "nonlinearity", "bad_pointing", "outliers",
+                             "angles_missing", "lu_eq_missing", "fresnel_angle_missing",
+                             "fresnel_default", "temp_variability_ed", "temp_variability_lu",
+                             "min_nbred", "min_nbrlu", "min_nbrlsky"]
+                    flagged = np.any([du.unpack_flags(datasetl1d['quality_flag'])[x] for x in flags], axis=0)
+                    ids = np.where(
+                        (datasetl1d['series_id'] == series_id[i]) & (flagged == False))
                     temp_arr[i] = np.mean(datasetl1d[variablestring].values[ids])
                 dataset_l2a[variablestring].values = temp_arr