docs: def _get_indexer(self, key):
          """ Get indexer for rasterio array.
  
          Parameter
          ---------
          key: tuple of int
  
          Returns
          -------
          band_key: an indexer for the 1st dimension
          window: two tuples. Each consists of (start, stop).
          squeeze_axis: axes to be squeezed
          np_ind: indexer for loaded numpy array
  
          See also
          --------
          indexing.decompose_indexer
          """
          assert len(key) == 3, 'rasterio datasets should always be 3D'
  
          # bands cannot be windowed but they can be listed
          band_key = key[0]
          np_inds = []
          # bands (axis=0) cannot be windowed but they can be listed
          if isinstance(band_key, slice):
              start, stop, step = band_key.indices(self.shape[0])
              band_key = np.arange(start, stop, step)
          # be sure we give out a list
          band_key = (np.asarray(band_key) + 1).tolist()
          if isinstance(band_key, list):  # if band_key is not a scalar
              np_inds.append(slice(None))
  
          # but other dims can only be windowed
          window = []
          squeeze_axis = []
          for i, (k, n) in enumerate(zip(key[1:], self.shape[1:])):
              if isinstance(k, slice):
                  # step is always positive. see indexing.decompose_indexer
                  start, stop, step = k.indices(n)
                  np_inds.append(slice(None, None, step))
              elif is_scalar(k):
                  # windowed operations will always return an array
                  # we will have to squeeze it later
                  squeeze_axis.append(- (2 - i))
                  start = k
                  stop = k + 1
              else:
                  start, stop = np.min(k), np.max(k) + 1
                  np_inds.append(k - start)
              window.append((start, stop))
  
          if isinstance(key[1], np.ndarray) and isinstance(key[2], np.ndarray):
              # do outer-style indexing
              np_inds[-2:] = np.ix_(*np_inds[-2:])
  
          return band_key, tuple(window), tuple(squeeze_axis), tuple(np_inds)