def isel_points(self, dim='points', **indexers):
        # type: (...) -> Dataset
        """Returns a new dataset with each array indexed pointwise along the
        specified dimension(s).

        This method selects pointwise values from each array and is akin to
        the NumPy indexing behavior of `arr[[0, 1], [0, 1]]`, except this
        method does not require knowing the order of each array's dimensions.

        Parameters
        ----------
        dim : str or DataArray or pandas.Index or other list-like object, optional
            Name of the dimension to concatenate along. If dim is provided as a
            string, it must be a new dimension name, in which case it is added
            along axis=0. If dim is provided as a DataArray or Index or
            list-like object, its name, which must not be present in the
            dataset, is used as the dimension to concatenate along and the
            values are added as a coordinate.
        **indexers : {dim: indexer, ...}
            Keyword arguments with names matching dimensions and values given
            by array-like objects. All indexers must be the same length and
            1 dimensional.

        Returns
        -------
        obj : Dataset
            A new Dataset with the same contents as this dataset, except each
            array and dimension is indexed by the appropriate indexers. With
            pointwise indexing, the new Dataset will always be a copy of the
            original.

        See Also
        --------
        Dataset.sel
        Dataset.isel
        Dataset.sel_points
        DataArray.isel_points
        """  # noqa
        warnings.warn('Dataset.isel_points is deprecated: use Dataset.isel()'
                      'instead.', DeprecationWarning, stacklevel=2)

        indexer_dims = set(indexers)

        def take(variable, slices):
            # Note: remove helper function when once when numpy
            # supports vindex https://github.com/numpy/numpy/pull/6075
            if hasattr(variable.data, 'vindex'):
                # Special case for dask backed arrays to use vectorised list
                # indexing
                sel = variable.data.vindex[slices]
            else:
                # Otherwise assume backend is numpy array with 'fancy' indexing
                sel = variable.data[slices]
            return sel

        def relevant_keys(mapping):
            return [k for k, v in mapping.items()
                    if any(d in indexer_dims for d in v.dims)]

        coords = relevant_keys(self.coords)
        indexers = [(k, np.asarray(v))  # type: ignore
                    for k, v in indexers.items()]
        indexers_dict = dict(indexers)
        non_indexed_dims = set(self.dims) - indexer_dims
        non_indexed_coords = set(self.coords) - set(coords)

        # All the indexers should be iterables
        # Check that indexers are valid dims, integers, and 1D
        for k, v in indexers:
            if k not in self.dims:
                raise ValueError("dimension %s does not exist" % k)
            if v.dtype.kind != 'i':  # type: ignore
                raise TypeError('Indexers must be integers')
            if v.ndim != 1:  # type: ignore
                raise ValueError('Indexers must be 1 dimensional')

        # all the indexers should have the same length
        lengths = set(len(v) for k, v in indexers)
        if len(lengths) > 1:
            raise ValueError('All indexers must be the same length')

        # Existing dimensions are not valid choices for the dim argument
        if isinstance(dim, str):
            if dim in self.dims:
                # dim is an invalid string
                raise ValueError('Existing dimension names are not valid '
                                 'choices for the dim argument in sel_points')

        elif hasattr(dim, 'dims'):
            # dim is a DataArray or Coordinate
            if dim.name in self.dims:
                # dim already exists
                raise ValueError('Existing dimensions are not valid choices '
                                 'for the dim argument in sel_points')

        # Set the new dim_name, and optionally the new dim coordinate
        # dim is either an array-like or a string
        if not utils.is_scalar(dim):
            # dim is array like get name or assign 'points', get as variable
            dim_name = 'points' if not hasattr(dim, 'name') else dim.name
            dim_coord = as_variable(dim, name=dim_name)
        else:
            # dim is a string
            dim_name = dim
            dim_coord = None  # type: ignore

        reordered = self.transpose(
            *(list(indexer_dims) + list(non_indexed_dims)))

        variables = OrderedDict()  # type: ignore

        for name, var in reordered.variables.items():
            if name in indexers_dict or any(
                    d in indexer_dims for d in var.dims):
                # slice if var is an indexer or depends on an indexed dim
                slc = [indexers_dict[k]
                       if k in indexers_dict
                       else slice(None) for k in var.dims]

                var_dims = [dim_name] + [d for d in var.dims
                                         if d in non_indexed_dims]
                selection = take(var, tuple(slc))
                var_subset = type(var)(var_dims, selection, var.attrs)
                variables[name] = var_subset
            else:
                # If not indexed just add it back to variables or coordinates
                variables[name] = var

        coord_names = (set(coords) & set(variables)) | non_indexed_coords

        dset = self._replace_vars_and_dims(variables, coord_names=coord_names)
        # Add the dim coord to the new dset. Must be done after creation
        # because_replace_vars_and_dims can only access existing coords,
        # not add new ones
        if dim_coord is not None:
            dset.coords[dim_name] = dim_coord
        return dset