fixing interpolate and other changes to make outdoor_test work and propagate the different types of uncertainties

Author: pdevis

examples/test_radiance_and_irradiance_libradtran.npy

@@ -67,6 +67,9 @@ def calibrate_l1a(self, measurandstring, dataset_l0, dataset_l0_bla, swir=False)
         if self.context.get_config_value("plot_l1a"):
             self.plot.plot_scans_in_series(measurandstring,dataset_l1a)
 
+        if self.context.get_config_value("plot_l1a_diff"):
+            self.plot.plot_diff_scans(measurandstring,dataset_l1a)
+
         return dataset_l1a
 
     def prepare_calibration_data(self,measurandstring, swir=False):

hypernets_processor/calibration/calibrate.py

@@ -4,6 +4,7 @@
 
 from hypernets_processor.version import __version__
 from hypernets_processor.data_utils.average import Average
+from hypernets_processor.plotting.plotting import Plotting
 from hypernets_processor.data_io.hypernets_writer import HypernetsWriter
 from hypernets_processor.data_io.data_templates import DataTemplates
 from hypernets_processor.combine_SWIR.measurement_functions.combine_factory import CombineFactory
@@ -28,6 +29,7 @@ def __init__(self,context,MCsteps=1000,parallel_cores=1):
         self.avg = Average(context=context)
         self.templ = DataTemplates(context)
         self.writer=HypernetsWriter(context)
+        self.plot=Plotting(context)
         self.context = context
 
     def combine(self,measurandstring,dataset_l1a,dataset_l1a_swir):
@@ -60,19 +62,19 @@ def combine(self,measurandstring,dataset_l1a,dataset_l1a_swir):
                      None]
         cov_systematic_input_qty_indep =  [None,
                      punpy.convert_corr_to_cov(dataset_l1b["corr_systematic_indep_" + measurandstring].values,
-                            dataset_l1b["u_systematic_indep_"+measurandstring].values),
+                            np.mean(dataset_l1b["u_systematic_indep_"+measurandstring].values,axis=1)),
                      None,
                      punpy.convert_corr_to_cov(dataset_l1b_swir["corr_systematic_indep_"+measurandstring].values,
-                            dataset_l1b_swir["u_systematic_indep_"+measurandstring].values),
+                            np.mean(dataset_l1b_swir["u_systematic_indep_"+measurandstring].values,axis=1)),
                      None]
         cov_systematic_input_qty_corr = [None,
                      punpy.convert_corr_to_cov(dataset_l1b["corr_systematic_corr_rad_irr_" + measurandstring].values,
-                            dataset_l1b["u_systematic_corr_rad_irr_"+measurandstring].values),
+                            np.mean(dataset_l1b["u_systematic_corr_rad_irr_"+measurandstring].values,axis=1)),
                      None,
                      punpy.convert_corr_to_cov(dataset_l1b_swir["corr_systematic_corr_rad_irr_"+measurandstring].values,
-                            dataset_l1b_swir["u_systematic_corr_rad_irr_"+measurandstring].values),
+                            np.mean(dataset_l1b_swir["u_systematic_corr_rad_irr_"+measurandstring].values,axis=1)),
                      None]
-
+        #todo do this more consistently with other modules, and do a direct copy for ranges that don't overlap
         dataset_l1b_comb = self.templ.l1b_template_from_combine(measurandstring,dataset_l1b,dataset_l1b_swir)
 
         self.process_measurement_function(measurandstring,dataset_l1b_comb,
@@ -84,15 +86,15 @@ def combine(self,measurandstring,dataset_l1a,dataset_l1a_swir):
                                           cov_systematic_input_qty_corr)
 
         if self.context.get_config_value("write_l1b"):
-            self.writer.write(dataset_l1b, overwrite=True)
+            self.writer.write(dataset_l1b_comb, overwrite=True)
 
         if self.context.get_config_value("plot_l1b"):
-            self.plot.plot_series_in_sequence(measurandstring,dataset_l1b)
+            self.plot.plot_series_in_sequence(measurandstring,dataset_l1b_comb)
 
         # if self.context.get_config_value("plot_diff"):
         #     self.plot.plot_diff_scans(measurandstring,dataset_l1a,dataset_l1b)
 
-        return dataset_l1b
+        return dataset_l1b_comb
 
     def find_input(self,variables,dataset):
         """
@@ -192,16 +194,16 @@ def process_measurement_function(self,measurandstring,dataset,measurement_functi
             warnings.simplefilter("ignore")
             u_random_measurand = self.prop.propagate_random(measurement_function,
                                                             input_quantities,
-                                                            u_random_input_quantities,
+                                                            u_random_input_quantities,param_fixed=[False,False,False,False,True],
                                                             repeat_dims=1)
             u_syst_measurand_indep,corr_syst_measurand_indep = self.prop.propagate_systematic(
                 measurement_function,input_quantities,
                 u_systematic_input_quantities_indep,
-                cov_x=cov_systematic_input_quantities_indep,return_corr=True,
+                cov_x=cov_systematic_input_quantities_indep,return_corr=True,param_fixed=[False,False,False,False,True],
                 repeat_dims=1,corr_axis=0)
             u_syst_measurand_corr,corr_syst_measurand_corr = self.prop.propagate_systematic(
                 measurement_function,input_quantities,u_systematic_input_quantities_corr,
-                cov_x=cov_systematic_input_quantities_corr,return_corr=True,
+                cov_x=cov_systematic_input_quantities_corr,return_corr=True,param_fixed=[False,False,False,False,True],
                 repeat_dims=1,corr_axis=0)
         dataset[measurandstring].values = measurand
         dataset["u_random_"+measurandstring].values = u_random_measurand

hypernets_processor/combine_SWIR/combine_SWIR.py

@@ -36,11 +36,9 @@ def average_l1b(self, measurandstring, dataset_l1a):
         dataset_l1b["corr_random_" + measurandstring].values = np.eye(
                 len(dataset_l1b["u_random_" + measurandstring].values))
         dataset_l1b["corr_systematic_indep_"+measurandstring].values = \
-            self.calc_mean_masked(dataset_l1a,"corr_systematic_indep_"+
-                                  measurandstring,corr=True)
+                dataset_l1a["corr_systematic_indep_"+measurandstring].values
         dataset_l1b["corr_systematic_corr_rad_irr_"+measurandstring].values = \
-            self.calc_mean_masked(dataset_l1a,"corr_systematic_corr_rad_irr_"+
-                                  measurandstring,corr=True)
+                dataset_l1a["corr_systematic_corr_rad_irr_"+measurandstring].values
 
         return dataset_l1b
 
@@ -49,6 +47,15 @@ def calc_mean_masked(self, dataset, var, rand_unc=False, corr=False):
         if corr:
             out = np.empty\
                 ((len(series_id), len(dataset['wavelength']), len(dataset['wavelength'])))
+
+            print("corr",dataset[var].values)
+            for i in range(len(series_id)):
+                ids = np.where((dataset['series_id'] == series_id[i]) & np.invert(
+                    DatasetUtil.unpack_flags(dataset["quality_flag"])["outliers"]))
+                out[i] = np.mean(dataset[var].values[:,:,ids],axis=3)[:,:,0]
+
+            out = np.mean(out, axis=0)
+
         else:
             out = np.empty((len(series_id), len(dataset['wavelength'])))
 
@@ -57,11 +64,8 @@ def calc_mean_masked(self, dataset, var, rand_unc=False, corr=False):
                            np.invert
                 (DatasetUtil.unpack_flags(dataset["quality_flag"])["outliers"]))
             out[i] = np.mean(dataset[var].values[:, ids], axis=2)[:, 0]
-            print("ids",i,ids,out[i].shape)
 
             if rand_unc:
                 out[i] = out[i] / len(ids[0])
-        if corr:
-            out = np.mean(out, axis=0)
-        print(i,out[0].shape,out.shape,out.T.shape)
+
         return out.T

hypernets_processor/data_utils/average.py

@@ -1,5 +1,5 @@
 [Processor]
-
+archive_directory: C:\Users\pdv\PycharmProjects\hypernets_processor\hypernets_processor\out\
 network: w
 
 [Wavelength_cal]

hypernets_processor/etc/processor.config

@@ -1,6 +1,8 @@
 [Processor]
 
 network: l
+archive_directory: C:\Users\pdv\PycharmProjects\hypernets_processor\hypernets_processor\out\
+version: 0.1
 
 [Wavelength_cal]
 mapping_vis_a: +1.6499700E+02
@@ -24,10 +26,10 @@ measurement_function_calibrate: StandardMeasurementFunction
 model=series_rep,series_id,vaa,azimuth_ref,vza,mode,action,it,scan_total,series_time
 
 [Interpolate]
-measurement_function_interpolate: WaterNetworkInterpolationLinear
+measurement_function_interpolate: LandNetworkInterpolationIrradianceLinear
 
 [SurfaceReflectance]
-measurement_function_surface_reflectance: WaterNetworkProtocol
+measurement_function_surface_reflectance: LandNetworkProtocol
 
 [CombineSWIR]
 combine_lim_wav: 1000
@@ -36,14 +38,15 @@ measurement_function_combine: StepCombine
 [Output]
 write_l0: False
 write_l1a: True
-write_l1b: False
-write_l1c: False
+write_l1b: True
+write_l1c: True
 write_l2a: True
 
 plot_l0: False
-plot_l1a: False
+plot_l1a: True
+plot_l1a_diff: False
 plot_l1b: True
-plot_l1c: False
+plot_l1c: True
 plot_l2a: True
 
 [Flags]

hypernets_processor/etc/processor_2.config

@@ -37,7 +37,7 @@ write_l2a: True
 plotting_format: png
 plot_l0: False
 plot_l1a: True
-plot_diff: True
+plot_l1a_diff: True
 plot_l1b: True
 plot_l1c: True
 plot_l2a: True

hypernets_processor/etc/processor_land_defaults.config

@@ -75,7 +75,7 @@ write_l2a: True
 plotting_format: png
 plot_l0: False
 plot_l1a: True
-plot_diff: True
+plot_l1a_diff: True
 plot_l1b: True
 plot_l1c: True
 plot_l2a: True

hypernets_processor/etc/processor_water_defaults.config

@@ -39,14 +39,15 @@ def interpolate_l1b_w(self,dataset_l1a_rad,dataset_l1b_uprad, dataset_l1b_irr):
         dataset_l1b["acquisition_time"] = dataset_l1a_rad["acquisition_time"].sel(scan=upscan)
         # is this correct????
         dataset_l1b["u_random_upwelling_radiance"] = dataset_l1a_rad["u_random_radiance"].sel(scan=upscan)
-        dataset_l1b["u_systematic_upwelling_radiance"] = dataset_l1a_rad["u_systematic_radiance"].sel(scan=upscan)
+        dataset_l1b["u_systematic_indep_upwelling_radiance"] = dataset_l1a_rad["u_systematic_indep_radiance"].sel(scan=upscan)
+        dataset_l1b["u_systematic_corr_rad_irr_upwelling_radiance"] = dataset_l1a_rad["u_systematic_corr_rad_irr_radiance"].sel(scan=upscan)
         dataset_l1b["corr_random_upwelling_radiance"] = dataset_l1a_rad["corr_random_radiance"]
-        dataset_l1b["corr_systematic_upwelling_radiance"] = dataset_l1a_rad["corr_systematic_radiance"]
+        dataset_l1b["corr_systematic_indep_upwelling_radiance"] = dataset_l1a_rad["corr_systematic_indep_radiance"]
+        dataset_l1b["corr_systematic_corr_rad_irr_upwelling_radiance"] = dataset_l1a_rad["corr_systematic_corr_rad_irr_radiance"]
 
-        print("interpolate sky radiance")
-        print(dataset_l1b)
+        self.context.logger.info("interpolate sky radiance")
         dataset_l1b=self.interpolate_skyradiance(dataset_l1b, dataset_l1b_uprad)
-        print("interpolate irradiances")
+        self.context.logger.info("interpolate irradiances")
         dataset_l1b=self.interpolate_irradiance(dataset_l1b, dataset_l1b_irr)
 
         if self.context.get_config_value("write_l1b"):
@@ -57,7 +58,6 @@ def interpolate_l1c(self,dataset_l1b_rad,dataset_l1b_irr):
 
 
         dataset_l1c=self.templ.l1c_from_l1b_dataset(dataset_l1b_rad)
-        print(dataset_l1b_rad["acquisition_time"])
         dataset_l1c["acquisition_time"].values = dataset_l1b_rad["acquisition_time"].values
 
         dataset_l1c["radiance"].values = dataset_l1b_rad["radiance"].values
@@ -83,12 +83,14 @@ def interpolate_irradiance(self,dataset_l1c,dataset_l1b_irr):
 
         acqui_irr = dataset_l1b_irr['acquisition_time'].values
         acqui_rad = dataset_l1c['acquisition_time'].values
+        print("times",acqui_rad,acqui_irr)
+        print(dataset_l1b_irr["u_systematic_indep_irradiance"].values.shape)
         cov_indep = punpy.convert_corr_to_cov(
             dataset_l1b_irr["corr_systematic_indep_irradiance"].values,
-            dataset_l1b_irr["u_systematic_indep_irradiance"].values)
+            np.mean(dataset_l1b_irr["u_systematic_indep_irradiance"].values,axis=1))
         cov_corr = punpy.convert_corr_to_cov(
             dataset_l1b_irr["corr_systematic_corr_rad_irr_irradiance"].values,
-            dataset_l1b_irr["u_systematic_corr_rad_irr_irradiance"].values)
+            np.mean(dataset_l1b_irr["u_systematic_corr_rad_irr_irradiance"].values,axis=1))
 
         dataset_l1c = self.process_measurement_function("irradiance",dataset_l1c,interpolation_function.function,
                                                         [acqui_rad,acqui_irr,dataset_l1b_irr['irradiance'].values],
@@ -106,10 +108,27 @@ def interpolate_skyradiance(self,dataset_l1c,dataset_l1a_skyrad):
         acqui_irr = dataset_l1a_skyrad['acquisition_time'].values
         acqui_rad = dataset_l1c['acquisition_time'].values
 
-        dataset_l1c = self.process_measurement_function("downwelling_radiance",dataset_l1c,interpolation_function.function,
-                                                        [acqui_rad,acqui_irr,dataset_l1a_skyrad['radiance'].values],
-                                                        [None,None,dataset_l1a_skyrad['u_random_radiance'].values],
-                                                        [None,None,dataset_l1a_skyrad['u_systematic_radiance'].values])
+        cov_indep = punpy.convert_corr_to_cov(
+            dataset_l1a_skyrad["corr_systematic_indep_radiance"].values,
+            np.mean(dataset_l1a_skyrad["u_systematic_indep_radiance"].values,axis=1))
+        cov_corr = punpy.convert_corr_to_cov(
+            dataset_l1a_skyrad["corr_systematic_corr_rad_irr_radiance"].values,
+            np.mean(dataset_l1a_skyrad["u_systematic_corr_rad_irr_radiance"].values,
+                    axis=1))
+
+        dataset_l1c = self.process_measurement_function("downwelling_radiance",dataset_l1c,
+                                                        interpolation_function.function,
+                                                        [acqui_rad,acqui_irr,
+                                                         dataset_l1a_skyrad[
+                                                             'radiance'].values],
+                                                        [None,None,dataset_l1a_skyrad[
+                                                            'u_random_radiance'].values],
+                                                        [None,None,dataset_l1a_skyrad[
+                                                            'u_systematic_indep_radiance'].values],
+                                                        [None,None,dataset_l1a_skyrad[
+                                                            'u_systematic_corr_rad_irr_radiance'].values],
+                                                        [None,None,cov_indep],
+                                                        [None,None,cov_corr])
         return dataset_l1c
 
     def process_measurement_function(self,measurandstring,dataset,measurement_function,
@@ -118,22 +137,49 @@ def process_measurement_function(self,measurandstring,dataset,measurement_functi
                                      u_systematic_input_quantities_corr,
                                      cov_systematic_input_quantities_indep,
                                      cov_systematic_input_quantities_corr):
+
+        # datashape = input_quantities[0].shape
+        # for i in range(len(input_quantities)):
+        #     if len(input_quantities[i].shape) > len(datashape):
+        #         datashape = input_quantities[0].shape
+        #
+        # for i in range(len(input_quantities)):
+        #     print(input_quantities[i].shape)
+        #     if len(input_quantities[i].shape) < len(datashape):
+        #         if input_quantities[i].shape[0]==datashape[1]:
+        #             input_quantities[i] = np.tile(input_quantities[i],(datashape[0],1))
+        #         else:
+        #             input_quantities[i] = np.tile(input_quantities[i],(datashape[1],1)).T
+        #     print(input_quantities[i].shape)
+        #
+        #     if u_random_input_quantities[i] is not None:
+        #         if len(u_random_input_quantities[i].shape) < len(datashape):
+        #             u_random_input_quantities[i] = np.tile(u_random_input_quantities[i], (datashape[1], 1)).T
+        #     if u_systematic_input_quantities_indep[i] is not None:
+        #         if len(u_systematic_input_quantities_indep[i].shape) < len(datashape):
+        #             u_systematic_input_quantities_indep[i] = np.tile(u_systematic_input_quantities_indep[i], (datashape[1], 1)).T
+        #     if u_systematic_input_quantities_corr[i] is not None:
+        #         if len(u_systematic_input_quantities_corr[i].shape) < len(datashape):
+        #             u_systematic_input_quantities_corr[i] = np.tile(u_systematic_input_quantities_corr[i], (datashape[1], 1)).T
+        print(*input_quantities)
         measurand = measurement_function(*input_quantities)
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             u_random_measurand = self.prop.propagate_random(measurement_function,
                                                             input_quantities,
                                                             u_random_input_quantities,
-                                                            repeat_dims=1)
+                                                        param_fixed=[False,True,True],
+                                                        repeat_dims=1)
+            print("shapes", cov_systematic_input_quantities_indep[2].shape,input_quantities[2].shape,u_systematic_input_quantities_indep[2].shape)
             u_syst_measurand_indep,corr_syst_measurand_indep = self.prop.propagate_systematic(
                 measurement_function,input_quantities,
                 u_systematic_input_quantities_indep,
-                cov_x=cov_systematic_input_quantities_indep,return_corr=True,
-                repeat_dims=1,corr_axis=0)
+                cov_x=['rand']*len(u_random_input_quantities),return_corr=True,
+                corr_axis=0,param_fixed=[False,True,True],repeat_dims=1)
             u_syst_measurand_corr,corr_syst_measurand_corr = self.prop.propagate_systematic(
                 measurement_function,input_quantities,u_systematic_input_quantities_corr,
-                cov_x=cov_systematic_input_quantities_corr,return_corr=True,
-                repeat_dims=1,corr_axis=0)
+                cov_x=['rand']*len(u_random_input_quantities),return_corr=True,
+                corr_axis=0,param_fixed=[False,True,True],repeat_dims=1)
         dataset[measurandstring].values = measurand
         dataset["u_random_"+measurandstring].values = u_random_measurand
         dataset["u_systematic_indep_"+measurandstring].values = u_syst_measurand_indep
@@ -146,9 +192,3 @@ def process_measurement_function(self,measurandstring,dataset,measurement_functi
             "corr_systematic_corr_rad_irr_"+measurandstring].values = corr_syst_measurand_corr
 
         return dataset
-
-
-
-
-
-