script_plot.py

# Below parameters can be changed
import uuid

import dask.array as da
import matplotlib.pyplot as plt
import numpy as np
import spatialdata_plot  # noqa: F401
from matplotlib.colors import Normalize
from spatialdata import bounding_box_query
from spatialdata.transformations import get_transformation

import harpy as hp
from harpy.image._image import _get_spatial_element

sdata = hp.datasets.resolve_example()


img_layer = "raw_image"
labels_layer = "segmentation_mask"


channels_to_plot = [
    0,
    # "CD14",
    # "CD3",
]

dpi = 200
figsize = (10, 10)
scale = "scale3"  # set to scale0 for full resolution, will be slower
crd = [0, 6000, 0, 5000]  # set to None if you do not want to crop


se = _get_spatial_element(sdata, layer=img_layer)
channels = se.c.data

assert set(channels_to_plot).issubset(channels), (
    f"Not all channels could be found in layer '{img_layer}'. Please choose from: {channels}."
)

img_layer_crop = f"{img_layer}_{uuid.uuid4()}"
labels_layer_crop = f"{labels_layer}_{uuid.uuid4()}"

if crd is not None:
    sdata[img_layer_crop] = bounding_box_query(
        sdata[img_layer],
        axes=["x", "y"],
        min_coordinate=[crd[0], crd[2]],
        max_coordinate=[crd[1], crd[3]],
        target_coordinate_system=list(get_transformation(se, get_all=True).keys())[0],
    )

    sdata[labels_layer_crop] = bounding_box_query(
        sdata[labels_layer],
        axes=["x", "y"],
        min_coordinate=[crd[0], crd[2]],
        max_coordinate=[crd[1], crd[3]],
        target_coordinate_system=list(get_transformation(se, get_all=True).keys())[0],
    )

if len(channels_to_plot) > 1:
    axes = plt.subplots(len(channels_to_plot), 1, figsize=figsize)[1].flatten()
else:
    axes = [plt.subplots(len(channels_to_plot), 1, figsize=figsize)[1]]

for c, ax in zip(channels_to_plot, axes, strict=False):
    c_id = np.where(channels == c)[0].item()
    vmax = da.percentile(
        se.data[c_id].flatten(), q=99
    ).compute()  # clip to 99% percentile
    norm = Normalize(vmax=vmax, clip=True)
    indices_channels_for_segmentation = [
        np.where(channels == _channel)[0].item() for _channel in channels_to_plot
    ]

    sdata.pl.render_images(
        img_layer if crd is None else img_layer_crop,
        channel=c,
        cmap="grey",
        norm=norm,
        scale=scale,
    ).pl.render_labels(
        labels_layer if crd is None else labels_layer_crop,
        scale=scale,
        fill_alpha=0.2,
        outline_alpha=0.3,
    ).pl.show(
        title=str(c),
        ax=ax,
        colorbar=False,
        dpi=dpi,
        return_ax=False,
    )
    ax.tick_params(
        left=False, bottom=False, labelleft=False, labelbottom=False
    )  # to remove ticks

if crd is not None:
    del sdata[img_layer_crop]
    del sdata[labels_layer_crop]