def colorize(self, colormap):
        """Colorize the current image using `colormap`.

        .. note::

            Works only on "L" or "LA" images.

        """

        if self.mode not in ("L", "LA"):
            raise ValueError("Image should be grayscale to colorize")

        if self.mode == "LA":
            alpha = self.data.sel(bands=['A'])
        else:
            alpha = None

        l_data = self.data.sel(bands=['L'])

        def _colorize(l_data, colormap):
            # 'l_data' is (1, rows, cols)
            # 'channels' will be a list of 3 (RGB) or 4 (RGBA) arrays
            channels = colormap.colorize(l_data)
            return np.concatenate(channels, axis=0)

        new_data = l_data.data.map_blocks(_colorize, colormap,
                                          chunks=(colormap.colors.shape[1],) + l_data.data.chunks[1:],
                                          dtype=np.float64)

        if colormap.colors.shape[1] == 4:
            mode = "RGBA"
        elif alpha is not None:
            new_data = da.concatenate([new_data, alpha.data], axis=0)
            mode = "RGBA"
        else:
            mode = "RGB"

        # copy the coordinates so we don't affect the original
        coords = dict(self.data.coords)
        coords['bands'] = list(mode)
        attrs = self.data.attrs
        dims = self.data.dims
        self.data = xr.DataArray(new_data, coords=coords, attrs=attrs, dims=dims)
Colorize the current image using `colormap`.

        .. note::

            Works only on "L" or "LA" images.