From aa1e76edb2c434c1a188952045c2078dc1a770bd Mon Sep 17 00:00:00 2001 From: ARamsay17 Date: Fri, 27 Oct 2023 11:21:28 +0100 Subject: [PATCH] ProcessInstrumentUncertainties.py, ProcessL2.py, RhoCorrection.py: L2 products now derived only for the wavelength subset used in rho calculation, uncertainties interpolated to same wavelength subset. --- Source/ProcessInstrumentUncertainties.py | 64 ++++++++++++------------ Source/ProcessL2.py | 12 ++--- Source/RhoCorrections.py | 3 ++ 3 files changed, 41 insertions(+), 38 deletions(-) diff --git a/Source/ProcessInstrumentUncertainties.py b/Source/ProcessInstrumentUncertainties.py index 90bdb91a..8410da75 100644 --- a/Source/ProcessInstrumentUncertainties.py +++ b/Source/ProcessInstrumentUncertainties.py @@ -487,7 +487,7 @@ def rrsHyperUNC(self, uncGrp: HDFGroup, rhoScalar: float, rhoVec: np.array, rhoD if rhoScalar is not None: # make rho a constant array if scalar rho = np.ones(len(list(esXstd.keys())))*rhoScalar - rhoUnc, _ = self.interp_common_wvls(np.array(rhoDelta, dtype=float), + rhoUNC, _ = self.interp_common_wvls(np.array(rhoDelta, dtype=float), waveSubset, np.asarray(list(esXstd.keys()), dtype=float)) else: @@ -545,7 +545,7 @@ def rrsHyperUNC(self, uncGrp: HDFGroup, rhoScalar: float, rhoVec: np.array, rhoD ones, ones] lw_uncertainties = [np.array(list(ltXstd.values())).flatten() * lt, - rhoUnc, + rhoUNC, np.array(list(liXstd.values())).flatten() * li, Cal['LI']/200, Cal['LT']/200, cStab['LI'], cStab['LT'], @@ -587,7 +587,7 @@ def rrsHyperUNC(self, uncGrp: HDFGroup, rhoScalar: float, rhoVec: np.array, rhoD output = {} # interpolate output uncertainties to the waveSubset (common wavebands of interpolated es, li, & lt) - wvls = np.asarray(list(xSlice['es'].keys()), dtype=float) + # wvls = np.asarray(list(xSlice['es'].keys()), dtype=float) lwAbsUnc, _ = self.interp_common_wvls(lwAbsUnc, np.array(uncGrp.getDataset("ES_RADCAL_CAL").columns['1'], dtype=float), @@ -604,34 +604,34 @@ def rrsHyperUNC(self, uncGrp: HDFGroup, rhoScalar: float, rhoVec: np.array, rhoD ## Band Convolution of Uncertainties if ConfigFile.settings["bL2WeightSentinel3A"]: - output["lwUNC_Sentinel3A"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_Sentinel3A"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "S3A") - output["rrsUNC_Sentinel3A"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_Sentinel3A"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "S3A") elif ConfigFile.settings["bL2WeightSentinel3B"]: - output["lwUNC_Sentinel3B"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_Sentinel3B"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "S3B") - output["rrsUNC_Sentinel3B"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_Sentinel3B"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "S3B") if ConfigFile.settings['bL2WeightMODISA']: - output["lwUNC_MODISA"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_MODISA"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "MOD-A") - output["rrsUNC_MODISA"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_MODISA"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "MOD-A") if ConfigFile.settings['bL2WeightMODIST']: - output["lwUNC_MODIST"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_MODIST"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "MOD-T") - output["rrsUNC_MODIST"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_MODIST"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "MOD-T") if ConfigFile.settings['bL2WeightVIIRSN']: - output["lwUNC_VIIRSN"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_VIIRSN"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "VIIRS") - output["rrsUNC_VIIRSN"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_VIIRSN"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "VIIRS") if ConfigFile.settings['bL2WeightVIIRSJ']: - output["lwUNC_VIIRSJ"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_VIIRSJ"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "VIIRS") - output["rrsUNC_VIIRSJ"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_VIIRSJ"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "VIIRS") pass output.update({"lwUNC": lwAbsUnc, "rrsUNC": rrsAbsUnc}) @@ -777,49 +777,49 @@ def rrsHyperUNCFACTORY(self, node, uncGrp, rhoScalar, rhoVec, rhoDelta, waveSubs output = {} # create dictionary to store uncertainty values which are returned from methods # interpolate output uncertainties to the waveSubset (common wavebands of interpolated es, li, & lt) - wvls = np.asarray(list(xSlice['es'].keys()), dtype=float) + # wvls = np.asarray(list(xSlice['es'].keys()), dtype=float) lwAbsUnc, _ = self.interp_common_wvls(lwAbsUnc, np.array(uncGrp.getDataset("ES_RADCAL_UNC").columns['wvl'], dtype=float), - wvls) + waveSubset) lw_vals, _ = self.interp_common_wvls(lw_vals, np.array(uncGrp.getDataset("ES_RADCAL_UNC").columns['wvl'], dtype=float), - wvls) + waveSubset) rrsAbsUnc, _ = self.interp_common_wvls(rrsAbsUnc, np.array(uncGrp.getDataset("ES_RADCAL_UNC").columns['wvl'], dtype=float), - wvls) + waveSubset) rrs_vals, _ = self.interp_common_wvls(rrs_vals, np.array(uncGrp.getDataset("ES_RADCAL_UNC").columns['wvl'], dtype=float), - wvls) + waveSubset) if ConfigFile.settings["bL2WeightSentinel3A"]: - output["lwUNC_Sentinel3A"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_Sentinel3A"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "S3A") - output["rrsUNC_Sentinel3A"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_Sentinel3A"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "S3A") elif ConfigFile.settings["bL2WeightSentinel3B"]: - output["lwUNC_Sentinel3B"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_Sentinel3B"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "S3B") - output["rrsUNC_Sentinel3B"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_Sentinel3B"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "S3B") if ConfigFile.settings['bL2WeightMODISA']: - output["lwUNC_MODISA"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_MODISA"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "MOD-A") - output["rrsUNC_MODISA"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_MODISA"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "MOD-A") if ConfigFile.settings['bL2WeightMODIST']: - output["lwUNC_MODIST"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_MODIST"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "MOD-T") - output["rrsUNC_MODIST"] = Convolve.band_Conv_Uncertainty([rrs_vals,wvls], + output["rrsUNC_MODIST"] = Convolve.band_Conv_Uncertainty([rrs_vals,waveSubset], [rrsAbsUnc, None], "MOD-T") if ConfigFile.settings['bL2WeightVIIRSN']: - output["lwUNC_VIIRSN"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_VIIRSN"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "VIIRS") - output["rrsUNC_VIIRSN"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_VIIRSN"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "VIIRS") if ConfigFile.settings['bL2WeightVIIRSJ']: - output["lwUNC_VIIRSJ"] = Convolve.band_Conv_Uncertainty([lw_vals, wvls], + output["lwUNC_VIIRSJ"] = Convolve.band_Conv_Uncertainty([lw_vals, waveSubset], [lwAbsUnc, None], "VIIRS") - output["rrsUNC_VIIRSJ"] = Convolve.band_Conv_Uncertainty([rrs_vals, wvls], + output["rrsUNC_VIIRSJ"] = Convolve.band_Conv_Uncertainty([rrs_vals, waveSubset], [rrsAbsUnc, None], "VIIRS") pass output.update({"lwUNC": lwAbsUnc, "rrsUNC": rrsAbsUnc}) diff --git a/Source/ProcessL2.py b/Source/ProcessL2.py index 895a3aae..76de22b5 100644 --- a/Source/ProcessL2.py +++ b/Source/ProcessL2.py @@ -373,9 +373,9 @@ def spectralReflectance(node, sensor, timeObj, xSlice, F0, F0_unc, rhoScalar, rh rrsUNC[k] = 0 deleteKey = [] - for k in esXSlice: # loop through wavebands - if (k in liXSlice) and (k in ltXSlice): - + for wvl in waveSubset: # loop through wavebands + k = str(wvl) + if (k in esXSlice) and (k in liXSlice) and (k in ltXSlice): # Initialize the new dataset if this is the first slice if k not in newESData.columns: newESData.columns[k] = [] @@ -440,7 +440,7 @@ def spectralReflectance(node, sensor, timeObj, xSlice, F0, F0_unc, rhoScalar, rh rrs = lw / es #Calculate the normalized water leaving radiance - nLw = rrs*f0 + nLw = rrs*f0 # need to chop the keys to match Z17 output # nLw uncertainty; nLwUNC[k] = np.power((rrsUNC[k]*f0)**2 + f0UNC**2, 0.5) @@ -460,7 +460,7 @@ def spectralReflectance(node, sensor, timeObj, xSlice, F0, F0_unc, rhoScalar, rh newLTDataMedian.columns[k].append(ltMedian) # Only populate valid wavelengths. Mark others for deletion - if float(k) in waveSubset: + if float(k) in waveSubset: # should be redundant! newRrsUncorrData.columns[k].append(rrs_uncorr) newLWData.columns[k].append(lw) newRrsData.columns[k].append(rrs) @@ -1532,7 +1532,7 @@ def _sliceRawData(ES_raw, LI_raw, LT_raw): if k != "Datetag" and k != "Datetime" and k != "Timetag2": wavelengthStr.append(k) wavelength.append(float(k)) - waveSubset = wavelength # Only used for Zhang; No subsetting for threeC or Mobley corrections + waveSubset = wavelength # Only used for Zhang; No subsetting for threeC or Mobley corrections rhoVec = {} Rho_Uncertainty_Obj = Propagate(M=100, cores=1) diff --git a/Source/RhoCorrections.py b/Source/RhoCorrections.py index 6190e1ee..751ad325 100644 --- a/Source/RhoCorrections.py +++ b/Source/RhoCorrections.py @@ -147,6 +147,9 @@ def ZhangCorr(windSpeedMean, AOD, cloud, sza, wTemp, sal, relAz, waveBands, Prop Utilities.writeLogFile(msg) # === environmental conditions during experiment === + if AOD > 0.2: # clip AOD to 0.2 to ensure no error in Z17, potential underestimation of uncertainty however + AOD = 0.2 + env = {'wind': windSpeedMean, 'od': AOD, 'C': cloud, 'zen_sun': sza, 'wtem': wTemp, 'sal': sal} # === The sensor ===