def _geoville_index_by_percentile(self, data, percentile):
        """ Calculate percentile of numpy stack and return the index of the chosen pixel. """
        # no_obs = bn.allnan(arr_tmp["data"], axis=0)
        data_tmp = np.array(data, copy=True)
        valid_obs = np.sum(np.isfinite(data_tmp), axis=0)
        # replace NaN with maximum
        max_val = np.nanmax(data_tmp) + 1
        data_tmp[np.isnan(data_tmp)] = max_val
        # sort - former NaNs will move to the end
        ind_tmp = np.argsort(data_tmp, kind="mergesort", axis=0)
        # desired position as well as floor and ceiling of it
        k_arr = (valid_obs - 1) * (percentile / 100.0)
        k_arr = np.where(k_arr < 0, 0, k_arr)
        f_arr = np.floor(k_arr + 0.5)
        f_arr = f_arr.astype(np.int)
        # get floor value of reference band and index band
        ind = f_arr.astype("int16")
        y_val, x_val = ind_tmp.shape[1], ind_tmp.shape[2]
        y_val, x_val = np.ogrid[0:y_val, 0:x_val]
        idx = np.where(valid_obs == 0, self.max_index, ind_tmp[ind, y_val, x_val])
        return idx