docs: 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