Merge remote-tracking branch 'origin/L1_and_L2_processing_structure' into L1_and_L2_processing_structure

# Conflicts:
#	hypernets_processor/data_io/format/variables.py
#	hypernets_processor/surface_reflectance/surface_reflectance.py

Author: pdevis

hypernets_processor/calibration/calibrate.py

@@ -19,7 +19,6 @@
 __email__ = "[email protected]"
 __status__ = "Development"
 
-
 class Calibrate:
     def __init__(self, context, MCsteps=1000, parallel_cores=0):
         self._measurement_function_factory = MeasurementFunctionFactory()
@@ -121,7 +120,7 @@ def calc_mean_masked(self, dataset, var, rand_unc=False, corr=False):
         else:
             out = np.empty((len(series_id), len(dataset['wavelength'])))
         for i in range(len(series_id)):
-            ids = np.where((dataset['series_id'] == series_id[i]) & (dataset['quality_flag'] == 1))
+            ids = np.where((dataset['series_id'] == series_id[i]) & (dataset['quality_flag'] == 0))
             out[i] = np.mean(dataset[var].values[:, ids], axis=2)[:, 0]
             if rand_unc:
                 out[i] = out[i] / len(ids[0])

hypernets_processor/data_io/hypernets_reader.py

@@ -132,7 +132,7 @@ def read_wavelength(self, pix):
         print("Wavelength range:", min(wvl), "-", max(wvl))
         return wvl
 
-    def read_series(self, seq_dir, series, lat, lon, metadata, fileformat):
+    def read_series(self, seq_dir, series, lat, lon, metadata, flag, fileformat):
 
         model_name = self.model
 
@@ -195,8 +195,7 @@ def read_series(self, seq_dir, series, lat, lon, metadata, fileformat):
 
         print("Wvl and Scan Dimensions:", len(wvl), scanDim)
         # use template from variables and metadata in format
-        ds = HypernetsDSBuilder()
-        ds = ds.create_ds_template(dim_sizes_dict=dim_sizes_dict, ds_format=fileformat)
+        ds = HypernetsDSBuilder().create_ds_template(dim_sizes_dict=dim_sizes_dict, ds_format=fileformat)
 
         ds["wavelength"] = wvl
         # ds["bandwidth"]=wvl
@@ -283,7 +282,7 @@ def read_series(self, seq_dir, series, lat, lon, metadata, fileformat):
                 ds["solar_zenith_angle"][scan_number] = get_altitude(float(lat), float(lon), acquisitionTime)
                 ds["solar_azimuth_angle"][scan_number] = get_azimuth(float(lat), float(lon), acquisitionTime)
 
-                # ds['quality_flag']
+                ds['quality_flag'][scan_number] = flag
                 ds['integration_time'][scan_number] = header['integration_time']
                 ds['temperature'][scan_number] = header['temperature']
 
@@ -317,7 +316,7 @@ def read_series(self, seq_dir, series, lat, lon, metadata, fileformat):
     def read_metadata(self, seq_dir):
 
         model_name=self.model
-
+        flag=0
         #     Spectra name : AA_BBB_CCCC_D_EEEE_FFF_GG_HHHH_II_JJJJ.spe
 
         #     A : iterator over "the sequence repeat time"
@@ -373,6 +372,7 @@ def read_metadata(self, seq_dir):
             if not lat.isdigit():
                 print("Latitude is not given, use default")
                 lat=self.context.get_config_value("lat")
+                flag=flag+2**self.context.get_config_value("lat_default")
 
             if 'longitude' in globalattr.keys():
                 lon = globalattr['longitude']
@@ -382,6 +382,7 @@ def read_metadata(self, seq_dir):
             if not lon.isdigit():
                 print("Longitude is not given, use default")
                 lon=self.context.get_config_value("lon")
+                flag=flag+2**self.context.get_config_value("lon_default")
 
             # 2. Estimate wavelengths - NEED TO CHANGE HERE!!!!!!
             # ----------------------
@@ -433,7 +434,7 @@ def read_metadata(self, seq_dir):
             print("Missing metadata file in sequence directory - check sequence directory")
             exit()
 
-        return seq, lat, lon, cc, metadata, seriesIrr, seriesRad, seriesBlack, seriesPict
+        return seq, lat, lon, cc, metadata, seriesIrr, seriesRad, seriesBlack, seriesPict,flag
 
     def read_sequence(self, seq_dir, setfile=None):
 
@@ -446,11 +447,11 @@ def read_sequence(self, seq_dir, setfile=None):
         L0_RAD = None
         L0_BLA = None
 
-        seq, lat, lon, cc, metadata, seriesIrr, seriesRad, seriesBlack, seriesPict = self.read_metadata(
+        seq, lat, lon, cc, metadata, seriesIrr, seriesRad, seriesBlack, seriesPict, flag = self.read_metadata(
             seq_dir)
 
         if seriesIrr:
-            L0_IRR = self.read_series(seq_dir, seriesIrr, lat, lon, metadata, "L0_IRR")
+            L0_IRR = self.read_series(seq_dir, seriesIrr, lat, lon, metadata, flag, "L0_IRR")
 
             # can't use this when non concatanted spectra
         #         if all([os.path.isfile(os.path.join(seq_dir,"RADIOMETER/",f)) for f in seriesIrr]):
@@ -461,7 +462,7 @@ def read_sequence(self, seq_dir, setfile=None):
             print("No irradiance data for this sequence")
 
         if seriesRad:
-            L0_RAD = self.read_series(seq_dir, seriesRad, lat, lon, metadata, "L0_RAD")
+            L0_RAD = self.read_series(seq_dir, seriesRad, lat, lon, metadata,flag, "L0_RAD")
         #         if all([os.path.isfile(os.path.join(seq_dir,"RADIOMETER/",f)) for f in seriesRad]):
         #             L0_RAD=read_series(seriesRad,cc, lat, lon, metadata, "L0_RAD")
         #         else:
@@ -470,7 +471,7 @@ def read_sequence(self, seq_dir, setfile=None):
             print("No radiance data for this sequence")
 
         if seriesBlack:
-            L0_BLA = self.read_series(seq_dir, seriesBlack, lat, lon, metadata, "L0_BLA")
+            L0_BLA = self.read_series(seq_dir, seriesBlack, lat, lon, metadata,flag, "L0_BLA")
             # if all([os.path.isfile(os.path.join(seq_dir, "RADIOMETER/", f)) for f in seriesBlack]):
             #     L0_BLA = self.read_series(seq_dir, seriesBlack, cc, lat, lon, metadata, "L0_BLA")
         else:
@@ -484,7 +485,6 @@ def read_sequence(self, seq_dir, setfile=None):
         return L0_IRR, L0_RAD, L0_BLA
 
     def read_flag(self, flagint):
-    #flag = 2 ** 1 + 2 ** 3 + 2 ** 5  # 2+8+32
         flagarray = 2 ** (np.linspace(0, 31, 32))
         flags = []
         while flagint:

hypernets_processor/etc/processor.config

@@ -57,3 +57,20 @@ test_sun_wave: 750
 test_sun_threshold: 0.05
 test_var_wave: 780
 test_var_threshold: 0.10
+
+[Flags]
+lon_default:7
+lat_default:8
+angles_missing:9
+lu_eq_missing:10
+fresnel_angle_missing:11
+fresnel_default:12
+temp_variability_ed:13
+temp_variability_lu:14
+## inf_or_nan:15
+min_nbred:16
+min_nbrlu:17
+min_nrblsky:18
+simil_fail:19
+wind_default:20
+

hypernets_processor/rhymer/rhymer/hypstar/rhymer_hypstar.py

@@ -85,22 +85,28 @@ def qc_scan(self, dataset, measurandstring):
 
                 ## continue if value exceeds the cv threshold
                 if v > self.context.get_config_value("diff_threshold"):
-                    # set quality flag to 10
-                    dataset['quality_flag'][scans[i]] = 10
+                    # get flag value for the temporal variability
+                    if measurandstring == 'irradiance':
+                        flagval=2 ** (self.context.get_config_value("temp_variability_ed"))
+                    else:
+                        flagval = 2 ** (self.context.get_config_value("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: print('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
-                    dataset['quality_flag'][scans[i]] = 10
-                    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']])))
+                # ## 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']])))
 
             return dataset
 
@@ -111,6 +117,9 @@ 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 = self.context.get_config_value("min_nbrlu")
+        flag_nbred = self.context.get_config_value("min_nbred")
+        flag_nbrlsky = self.context.get_config_value("min_nbrlsky")
 
         if protocol != 'water_std':
             # here we should simply provide surface reflectance?
@@ -131,6 +140,8 @@ 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: print(
                         'NULL angles: Aquisition time {}, {}'.format(ts, ', '.join(
                             ['{}:{}'.format(k, scani[k].values) for k in ['scan', 'quality_flag']])))
@@ -155,7 +166,8 @@ 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:
-                    lu["quality_flag"][lu["viewing_azimuth_angle"] == i] = 10
+                    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
                     if self.context.get_config_value("verbosity") > 2:
                         ts = [datetime.utcfromtimestamp(x) for x in
                               lu['acquisition_time'][lu["viewing_azimuth_angle"] == i].values]
@@ -169,7 +181,8 @@ def cycleparse(self, rad, irr):
             senz_lsky = 180 - np.unique(lsky["viewing_zenith_angle"].values)
             for i in senz_lu:
                 if i not in senz_lsky:
-                    lu["quality_flag"][lu["viewing_zenith_angle"] == i] = 10
+                    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]
@@ -180,14 +193,21 @@ def cycleparse(self, rad, irr):
                                      in ['scan', 'quality_flag']])))
 
             # check if correct number of radiance and irradiance data
-            if lu.scan[lu['quality_flag'] <= 1].count() < nbrlu:
+            flag_nbrlu = self.context.get_config_value("min_nbrlu")
+            flag_nbred = self.context.get_config_value("min_nbred")
+            flag_nbrlsky = self.context.get_config_value("min_nbrlsky")
+
+            if lu.scan[lu['quality_flag'] <= 0].count() < nbrlu:
+                lu["quality_flag"].values = [lu.sel(scan=i)["quality_flag"]+2**flag_nbrlu for i in lu['scan']]
                 if self.context.get_config_value("verbosity") > 2:
                     print("No enough upwelling radiance data for sequence {}".format(lu.attrs['sequence_id']))
             if lsky.scan[lsky['quality_flag'] <= 1].count() < nbrlsky:
+                lsky["quality_flag"] = [lsky.sel(scan=i)["quality_flag"]+2**flag_nbrlsky for i in lsky['scan']]
                 if self.context.get_config_value("verbosity") > 2:
                     print("No enough downwelling radiance data for sequence {}".format(lsky.attrs['sequence_id']))
 
             if irr.scan[irr['quality_flag'] <= 1].count() < nbred:
+                irr["quality_flag"] = [irr.sel(scan=i)["quality_flag"] + 2 ** flag_nbred for i in irr['scan']]
                 if self.context.get_config_value("verbosity") > 2:
                     print("No enough irradiance data for sequence {}".format(irr.attrs['sequence_id']))
 
@@ -198,9 +218,11 @@ def get_wind(self, l1b):
         lat = l1b.attrs['site_latitude']
         lon = l1b.attrs['site_latitude']
         wind = []
+        flagval = self.context.get_config_value("wind_default")
         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
                 print("Default wind speed {}".format(self.context.get_config_value("wind_default")))
                 wind.append(self.context.get_config_value("wind_default"))
             else:
@@ -237,15 +259,15 @@ def get_fresnelrefl(self, l1b):
 
             ## get fresnel reflectance
             if self.context.get_config_value("fresnel_option") == 'Mobley':
-                # check here if no inconsistency between LUT and fresnel option??
-                rholut = self.rhymerproc.mobley_lut_read()
                 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])
                 else:
                     # add a quality flag!
-                    fresnel = 'fixed'
+                    flagval = self.context.get_config_value("fresnel_default")
+                    l1b["quality_flag"][i].values = l1b["quality_flag"][i].values + 2 ** flagval
+                    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")
@@ -262,31 +284,31 @@ def get_fresnelrefl(self, l1b):
 
         return l1b
 
-    def get_epsilon(self, rhow_nosc, wavelength):
-
-        # wavelength = l1b['wavelength'].values
-        epsilon = np.zeros(len(rhow_nosc))
-
-        ## compute similarity epsilon
-        for i in range(len(rhow_nosc)):
-            fail_simil, eps = self.qc_similarity(wavelength, rhow_nosc[i],
-                                                 self.similarity_wr,
-                                                 self.similarity_wp,
-                                                 self.similarity_w1,
-                                                 self.similarity_w2,
-                                                 self.similarity_alpha)
-
-            ## R2005 quality control
-            ## skip spectra not following similarity
-            if self.similarity_test:
-                if fail_simil:
-                    if verbosity > 2: print('Failed simil test.')
-                    continue
-                else:
-                    if verbosity > 2: print('Passed simil test.')
-            epsilon[i] = eps
-        # l1b["epsilon"].values=epsilon
-        return epsilon
+    # def get_epsilon(self, rhow_nosc, wavelength):
+    #
+    #     # wavelength = l1b['wavelength'].values
+    #     epsilon = np.zeros(len(rhow_nosc))
+    #
+    #     ## compute similarity epsilon
+    #     for i in range(len(rhow_nosc)):
+    #         fail_simil, eps = self.qc_similarity(wavelength, rhow_nosc[i],
+    #                                              self.similarity_wr,
+    #                                              self.similarity_wp,
+    #                                              self.similarity_w1,
+    #                                              self.similarity_w2,
+    #                                              self.similarity_alpha)
+    #
+    #         ## R2005 quality control
+    #         ## skip spectra not following similarity
+    #         if self.similarity_test:
+    #             if fail_simil:
+    #                 if verbosity > 2: print('Failed simil test.')
+    #                 continue
+    #             else:
+    #                 if verbosity > 2: print('Passed simil test.')
+    #         epsilon[i] = eps
+    #     # l1b["epsilon"].values=epsilon
+    #     return epsilon
 
     def get_rhow_nosc(self, l1b):
 
@@ -407,16 +429,16 @@ def fresnelrefl_qc_simil(self, l1b, wind):
                                                  self.similarity_w2,
                                                  self.similarity_alpha)
 
-            ## R2005 quality control
-            ## skip spectra not following similarity
-            if self.similarity_test:
-                if fail_simil:
-                    if verbosity > 2: print('Failed simil test.')
-                    simil_flag[i] = 10
-                    continue
-                else:
-                    if verbosity > 2: print('Passed simil test.')
-                    simil_flag[i] = 0
+            # ## R2005 quality control
+            # ## skip spectra not following similarity
+            # if self.similarity_test:
+            #     if fail_simil:
+            #         if verbosity > 2: print('Failed simil test.')
+            #         simil_flag[i] = 10
+            #         continue
+            #     else:
+            #         if verbosity > 2: print('Passed simil test.')
+            #         simil_flag[i] = 0
 
             ## R2005 correction
             if self.similarity_correct:

hypernets_processor/surface_reflectance/measurement_functions/water_network_protocol.py

@@ -13,7 +13,7 @@ def function(self,upwelling_radiance,downwelling_radiance,irradiance,rhof,epsilo
 
         rhow_all = [r - epsilon for r in rhow_nosc_all]
 
-        return lw_all, rhow_nosc_all, rhow_all
+        return lw_all, rhow_nosc_all, rhow_all, epsilon
 
     @staticmethod
     def get_name():