Perform histogram matching of the time-series with respect to a reference scene

            :param eopatch: eopatch holding the time-series and reference data
            :type eopatch: EOPatch
            :return: The same eopatch instance with the normalised time-series

def execute(self, eopatch):
        """ Perform histogram matching of the time-series with respect to a reference scene

            :param eopatch: eopatch holding the time-series and reference data
            :type eopatch: EOPatch
            :return: The same eopatch instance with the normalised time-series
        """
        feature_type, feature_name, new_feature_name = next(self.feature(eopatch))
        reference_type, reference_name = next(self.reference(eopatch))

        reference_scene = eopatch[reference_type][reference_name]
        # check if band dimension matches
        if reference_scene.shape[-1] != eopatch[feature_type][feature_name].shape[-1]:
            raise ValueError('Time-series and reference scene must have corresponding bands')

        eopatch[feature_type][new_feature_name] = np.zeros_like(eopatch[feature_type][feature_name])
        for source_id, source in enumerate(eopatch[feature_type][feature_name]):
            # mask-out same invalid pixels
            src_masked = np.where(np.isnan(reference_scene), np.nan, source)
            ref_masked = np.where(np.isnan(source), np.nan, reference_scene)
            # compute statistics
            std_ref = np.nanstd(ref_masked, axis=(0, 1), dtype=np.float64)
            std_src = np.nanstd(src_masked, axis=(0, 1), dtype=np.float64)
            mean_ref = np.nanmean(ref_masked, axis=(0, 1), dtype=np.float64)
            mean_src = np.nanmean(src_masked, axis=(0, 1), dtype=np.float64)
            # normalise values
            eopatch[feature_type][new_feature_name][source_id] = \
                source * (std_ref / std_src) + (mean_ref - (mean_src * (std_ref / std_src)))

        return eopatch