Failing to recreate Xenium Exlorer-style multi-stain segmentation panel plot

Hi, I am attempting to use Spatiadata's pl.render_images to recreate a Xenium Explorer-style plot for Xenium with the multi-stain segmentation panel. See the example here:

I've processed the data so far as:

import spatialdata as sd
import spatialdata_plot
from spatialdata_io import xenium
from spatialdata.models import get_channel_names
import scanpy as sc
from pathlib import Path
import shutil
import os
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.patches as patches
from matplotlib import pyplot as plt
path = Path().resolve()

# Set root directory
root_dir = ... # Skipping
# Set figure output directory
figdir = os.path.join(
    root_dir,
    'figures'
)
# Set directory to raw Xenium output bundles
raw_data_dir = os.path.join(
    root_dir,
    'data/raw/xenium'
)
# Set path to reference NicheCompass Anndata file
adata_file = os.path.join(
    root_dir,
    'output/nichecompass',
    'FB_SPI_mouse_filt_cells_NicheCompass_annotated__30012026.h5ad'
)
# Mapping Xenium slide names to Xenium output directories
xenium_output_dir = {
    '7dpi-C': 'output-XETG00155__0043039__7dpi-C__20240725__115033',
    'Sham-A': 'output-XETG00155__0043039__Sham-A__20240725__115033',
    '1dpi-C': 'output-XETG00155__0043039__1dpi-C__20240725__115033',
    '3dpi-A': 'output-XETG00155__0043042__3dpi-A__20240725__115033',
    '7dpi-B': 'output-XETG00335__0027406__7dpi-B__20240725__114550',
    '14dpi-B': 'output-XETG00335__0027281__14dpi-B__20240725__114550',
    '28dpi-A': 'output-XETG00155__0043039__28dpi-A__20240725__115033',
}
slide_xenium_data_dir = {}
for slide in xenium_output_dir.keys():
    slide_xenium_data_dir[slide] = os.path.join(
        raw_data_dir,
        xenium_output_dir[slide]
    )

# Load NicheCompass results
adata_nc = sc.read_h5ad(adata_file)

# Load Xenium data as SpatialData
slide = '7dpi-C'
sdata = xenium(
    path=slide_xenium_data_dir[slide],
    n_jobs=8,
    cell_boundaries=True,
    nucleus_boundaries=False,
    morphology_focus=True,
    cells_as_circles=False,
)

# Subset NicheCompass to slide
adata_nc_subset = adata_nc[adata_nc.obs['regionName']==slide,:].copy()

# Fix 8-digit hex-codes to 6-digit
adata_nc_subset.uns['Niches_colors'] = np.array(
    [c[:-2] for c in adata_nc_subset.uns['Niches_colors']]
)
# Copy over to Spatialdata object
sdata.tables['table'].obs_names = sdata.tables['table'].obs['cell_id']
sdata.tables['table'].obs['Niches'] = adata_nc_subset.obs['Niches']
sdata.tables['table'].obs['cell_type_coarse'] = adata_nc_subset.obs['cell_type_coarse']
sdata.tables['table'].obs['cell_type_granular'] = adata_nc_subset.obs['cell_type_granular']
sdata.tables['table'].uns['Niches_colors'] = adata_nc_subset.uns['Niches_colors']
sdata.tables['table'].uns['cell_type_coarse_colors'] = adata_nc_subset.uns['cell_type_coarse_colors']
sdata.tables['table'].uns['cell_type_granular_colors'] = adata_nc_subset.uns['cell_type_granular_colors']

# Reset index of .obs
sdata.tables['table'].obs.reset_index(drop=True, inplace=True)

# Reset table annotation to 'cell_boundaries'
sdata["table"].obs["region"] = "cell_boundaries"
sdata.set_table_annotates_spatialelement("table", region="cell_boundaries")

# Replace NA values
sdata.tables['table'].obs['Niches'] = (
    sdata.tables['table'].obs['Niches'].cat.add_categories(['Unknown'])
          .fillna('Unknown')
)

# Process colour palette
niches_g = sdata.tables['table'].obs['Niches'].cat.categories.tolist()
niches_p = sdata.tables['table'].uns['Niches_colors'].tolist()
# Add a separate value for cells lacking 'Niche' annotation
niches_p.append('#ffffff')

INFO     reading                                                                                                   
         [/lustre/scratch124/cellgen/bayraktar/kr23/projects/FB_SPI/data/raw/xenium/output-XETG00155__0043039__7dpi-](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/lustre/scratch124/cellgen/bayraktar/kr23/projects/FB_SPI/data/raw/xenium/output-XETG00155__0043039__7dpi-)
         C__20240725__115033/cell_feature_matrix.h5                                                                
<timed exec>:3: DeprecationWarning: `cell_boundaries` is being deprecated as an argument to `xenium.xenium` in SpatialData version 0.1, switch to `cells_boundaries` instead.
[/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata_io/_utils.py:48](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata_io/_utils.py:48): UserWarning: The "cell_id" column in the cells metadata table does not match the "cell_id" column in the annotation table. This could be due to trying to read a new version that is not supported yet. Please report this issue.
  return f(*args, **kwargs)
[/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata_io/_utils.py:48](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata_io/_utils.py:48): UserWarning: The cell_id column in the cell_labels_table does not match the cell_id column derived from the cell labels data. This could be due to trying to read a new version that is not supported yet. Please report this issue.
  return f(*args, **kwargs)
CPU times: user 23.3 s, sys: 4.31 s, total: 27.6 s
Wall time: 25.1 s
[/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata/_core/spatialdata.py:511](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata/_core/spatialdata.py:511): UserWarning: Converting `region_key: region` to categorical dtype.
  convert_region_column_to_categorical(table)

# Bounding box to region of interest
bb_xmin = 11500
bb_ymin = 4000
bb_w = 4000
bb_h = 4000
bb_xmax = bb_xmin + bb_w
bb_ymax = bb_ymin + bb_h

# Highlight region of interest
f, ax = plt.subplots(figsize=(12,24))
sdata.pl.render_shapes(
    'cell_boundaries',
    method = 'matplotlib',
    color = 'Niches',
    groups = niches_g,
    palette = niches_p,
    alpha_fill = 1
    ).pl.show(
        ax = ax,
        title = f'{slide} (zoom)',
    )
rect = patches.Rectangle((bb_xmin, bb_ymin), bb_w, bb_h, linewidth=1, edgecolor="red", facecolor="none")
ax.add_patch(rect)

# Query the bounding box
sdata_xenium_subset = sdata.query.bounding_box(
    axes=["x", "y"],
    min_coordinate=[bb_xmin, bb_ymin],
    max_coordinate=[bb_xmax, bb_ymax],
    target_coordinate_system="global",
)

sdata_xenium_subset.pl.render_shapes(
    'cell_boundaries',
    method = 'matplotlib',
    color = 'Niches',
    groups = niches_g,
    palette = niches_p,
    fill_alpha = 0.8
    ).pl.show(
        title = f'{slide} (zoom) (niches)',
        figsize = (12,12)
    )

The full tissue slide:

And the queried bounding box:

Now I am aware of this issue: https://github.com/scverse/spatialdata-plot/issues/370,
but don't even seem to get a similar result:

sdata_xenium_subset.pl.render_images(
    "morphology_focus",
    palette = ["#0F73E6", "#F300A5", "#A4A400", "#008A00"],
    channel = ['DAPI', 'ATP1A1/CD45/E-Cadherin', '18S', 'AlphaSMA/Vimentin']
).pl.show(
    title="Morphology image", 
    coordinate_systems="global",
    colorbar = False,
    figsize = (6,6)
)

[/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/tifffile/tifffile.py:8797](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/tifffile/tifffile.py:8797): UserWarning: <tifffile.TiffPage 0 @16> reading array from closed file
  warnings.warn(
INFO     Your image has 4 channels. Sampling categorical colors and using multichannel strategy 'stack' to render.

Individual channels look fine:

sdata_xenium_subset.pl.render_images(
    'morphology_focus',
    palette = ['#0f73e6'],
    channel = ['DAPI'],
    alpha = 1
    # norm = matplotlib.colors.Normalize(0, 5000, clip=True)
    ).pl.show(
        title = f'{slide} (zoom) (DAPI)',
        colorbar = False,
        figsize = (6,6)
    )

sdata_xenium_subset.pl.render_images(
    'morphology_focus',
    palette = ['#a4a400'],
    channel = ['18S'],
    # norm = matplotlib.colors.Normalize(0, 5000, clip=True)
    ).pl.show(
        title = f'{slide} (zoom) (18S)',
        colorbar = False,
        figsize = (6,6)
    )

sdata_xenium_subset.pl.render_images(
    'morphology_focus',
    palette = ['#f300a5'],
    channel = ['ATP1A1/CD45/E-Cadherin'],
    # norm = matplotlib.colors.Normalize(0, 8000, clip=True)
    ).pl.show(
        title = f'{slide} (zoom) (ATP1A1/CD45/E-Cadherin)',
        colorbar = False,
        figsize = (6,6)
    )

However none of the options to do normalisation clipping seem to affect any of the combined plots, neither when I do this for instance:

###
sdata_xenium_subset.pl.render_images(
    'morphology_focus',
    palette = [
        '#0f73e6',
        '#a4a400',
        '#f300a5'
        ],
    channel = [
        'DAPI',
        '18S',
        'ATP1A1/CD45/E-Cadherin'
        ],
    # norm = matplotlib.colors.Normalize(0, 8000, clip=True)
    ).pl.show(
        title = f'{slide} (zoom) (DAPI - 18S - E-Cadherin)',
        colorbar = False,
        figsize = (6,6),
    )

Is exactly the same look regardless of setting e.g. norm = matplotlib.colors.Normalize(0, 10000, clip=True), norm = matplotlib.colors.Normalize(0, 1.5, clip=True), or norm = matplotlib.colors.Normalize(0, 2000, clip=True).

This is what is found in the queried object:

sdata_xenium_subset

SpatialData object
├── Images
│     └── 'morphology_focus': DataTree[cyx] (4, 4000, 4000), (4, 2000, 2000), (4, 1000, 1000), (4, 500, 500), (4, 250, 250)
├── Labels
│     ├── 'cell_labels': DataTree[yx] (4000, 4000), (2000, 2000), (1000, 1000), (500, 500), (250, 250)
│     └── 'nucleus_labels': DataTree[yx] (4000, 4000), (2000, 2000), (1000, 1000), (500, 500), (250, 250)
├── Points
│     └── 'transcripts': DataFrame with shape: (<Delayed>, 13) (3D points)
├── Shapes
│     └── 'cell_boundaries': GeoDataFrame shape: (4500, 1) (2D shapes)
└── Tables
      └── 'table': AnnData (4500, 474)
with coordinate systems:
    ▸ 'global', with elements:
        morphology_focus (Images), cell_labels (Labels), nucleus_labels (Labels), transcripts (Points), cell_boundaries (Shapes)

and images:

sdata_xenium_subset.images
{'morphology_focus': <xarray.DataTree>
Group: [/](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/)
├── Group: [/scale0](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/scale0)
│       Dimensions:  (c: 4, y: 4000, x: 4000)
│       Coordinates:
│         * c        (c) <U22 352B 'DAPI' ... 'AlphaSMA/Vimentin'
│         * y        (y) float64 32kB 0.5 1.5 2.5 3.5 ... 3.998e+03 3.998e+03 4e+03
│         * x        (x) float64 32kB 0.5 1.5 2.5 3.5 ... 3.998e+03 3.998e+03 4e+03
│       Data variables:
│           image    (c, y, x) uint16 128MB dask.array<chunksize=(2, 4000, 4000), meta=np.ndarray>
├── Group: [/scale1](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/scale1)
│       Dimensions:  (c: 4, y: 2000, x: 2000)
│       Coordinates:
│         * c        (c) <U22 352B 'DAPI' ... 'AlphaSMA/Vimentin'
│         * y        (y) float64 16kB 1.0 3.0 5.0 7.0 ... 3.995e+03 3.997e+03 3.999e+03
│         * x        (x) float64 16kB 1.0 3.0 5.0 7.0 ... 3.995e+03 3.997e+03 3.999e+03
│       Data variables:
│           image    (c, y, x) uint16 32MB dask.array<chunksize=(2, 2000, 2000), meta=np.ndarray>
├── Group: [/scale2](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/scale2)
│       Dimensions:  (c: 4, y: 1000, x: 1000)
│       Coordinates:
│         * c        (c) <U22 352B 'DAPI' ... 'AlphaSMA/Vimentin'
│         * y        (y) float64 8kB 2.0 6.0 10.0 14.0 ... 3.99e+03 3.994e+03 3.998e+03
│         * x        (x) float64 8kB 2.0 6.0 10.0 14.0 ... 3.99e+03 3.994e+03 3.998e+03
│       Data variables:
...
          * y        (y) float64 2kB 8.0 24.0 40.0 56.0 ... 3.96e+03 3.976e+03 3.992e+03
          * x        (x) float64 2kB 8.0 24.0 40.0 56.0 ... 3.96e+03 3.976e+03 3.992e+03
        Data variables:
            image    (c, y, x) uint16 500kB dask.array<chunksize=(4, 250, 250), meta=np.ndarray>}
Output is truncated. View as a [scrollable element](command:cellOutput.enableScrolling?a43a3130-bed1-4d79-94bd-0ac7729b722a) or open in a [text editor](command:workbench.action.openLargeOutput?a43a3130-bed1-4d79-94bd-0ac7729b722a). Adjust cell output [settings](command:workbench.action.openSettings?%5B%22%40tag%3AnotebookOutputLayout%22%5D)...

Are you familiar with any ways to improve this plotting? I would love to have the custom control that Spatialdata offers to generate high-quality figures for publications, but so far the different scales of stains doesn't seem to work well.

Finally, I have also tried this issue: https://github.com/scverse/spatialdata-plot/issues/460, but got exactly the same issue in the end as the author mentioned in their last comment.

All my code was run with:

print(sd.__version__)
print(spatialdata_io.__version__)
print(spatialdata_plot.__version__)
0.5.0
0.3.0
0.2.12

I hope you can help me out!

Best,
Koen

Empieza reproduciendo el problema con las llamadas a render_images sobre la imagen morphology_focus consultada, comparando los canales individuales con el renderizado multicanal de stack. Lee cómo se aplican los argumentos palette, channel y norm en la ruta multicanal. Se considerará terminado cuando los canales Xenium combinados respondan a la normalización y se rendericen con los colores y el contraste solicitados.