Plotting#

This example shows how to load an MDHistoWorkspace, integrate a thin slice, extract the plotted data, build a projection-aware transform, and generate a publication-style 2D plot using Matplotlib and axisartist.

The helper functions are intended for normalized reciprocal-space data loaded with Mantid’s LoadMD and plotted after integrating one reciprocal-space axis over a finite thickness.

Overview#

The workflow is:

Load the MD workspace if it is not already present.
Integrate a thin slice along one axis.
Extract the 2D coordinates, signal, and errors from the sliced workspace.
Construct a projection transform from the oriented lattice and W matrix.
Plot the slice with transformed axes and percentile-based color scaling.

Example script#

import os

from mantid.simpleapi import LoadMD, IntegrateMDHistoWorkspace, mtd

import numpy as np
import scipy.linalg

import matplotlib.pyplot as plt
from matplotlib.transforms import Affine2D

from mpl_toolkits.axisartist import Axes, GridHelperCurveLinear
from mpl_toolkits.axisartist.grid_finder import (
    ExtremeFinderSimple,
    MaxNLocator,
)

# workspace
folder = "/SNS/EXAMPLES/TOPAZ/IPTS-36169/shared/nxv/"
filename = "Si_AG_300K_normalization/Si_AG_300K_(h,k,0)_[0,0,l]_[-10.0,10.0]_[-10.0,10.0]_[-10.0,10.0]_201x201x201_m-3m_sub_bkg.nxs"
ws = "Si_300K"

# slice parameters
value = 0.0
thickness = 0.1
normal = [0, 0, 1]

# colorbar limits
pmin = 0
pmax = 95


def load_data(filename, ws):

    if not mtd.doesExist(ws):
        LoadMD(
            Filename=os.path.join(folder, filename),
            LoadHistory=False,
            OutputWorkspace=ws,
        )


def make_slice(ws, thickness, value, normal=[0, 0, 1]):

    integration = [value - thickness, value + thickness]

    IntegrateMDHistoWorkspace(
        InputWorkspace=ws,
        P1Bin=integration if normal[0] == 1 else None,
        P2Bin=integration if normal[1] == 1 else None,
        P3Bin=integration if normal[2] == 1 else None,
        OutputWorkspace=ws + "_slice",
    )

    return integration


def extract_data(ws):

    dims = mtd[ws].getNonIntegratedDimensions()

    xs = np.meshgrid(
        *[
            np.linspace(
                dim.getMinimum() + dim.getBinWidth() / 2,
                dim.getMaximum() - dim.getBinWidth() / 2,
                dim.getNBins(),
            )
            for dim in dims
        ],
        indexing="ij"
    )

    signal = mtd[ws].getSignalArray().squeeze()
    errors = np.sqrt(mtd[ws].getErrorSquaredArray().squeeze())

    return xs, signal, errors


def projection_transform(ws, y, normal=[0, 0, 1]):

    ei = mtd[ws].getExperimentInfo(0)

    dims = mtd[ws].getNonIntegratedDimensions()

    names = np.array([dim.name for dim in dims])

    B = ei.sample().getOrientedLattice().getB().copy()
    W = ei.run().getLogData("W_MATRIX").value.reshape(3, 3)

    Bp = np.dot(B, W)

    Q, R = scipy.linalg.qr(Bp)

    ind = np.abs(normal) != 1

    titles = names[ind].tolist() + names[~ind].tolist()

    v = scipy.linalg.cholesky(np.dot(R.T, R)[ind][:, ind], lower=False)

    v /= v[0, 0]

    T = np.eye(3)
    T[:2, :2] = v

    s = np.diag(T).copy()
    T[1, 1] = 1

    T[0, 2] = -T[0, 1] * y.min()

    return T, s[1], titles


def make_plot(x, y, z, T, aspect, titles, integration, pmin=5, pmax=95):

    transform = Affine2D(T)

    extreme_finder = ExtremeFinderSimple(20, 20)

    grid_locator1 = MaxNLocator(nbins=10)
    grid_locator2 = MaxNLocator(nbins=10)

    grid_locator1.set_params(integer=True)
    grid_locator2.set_params(integer=True)

    minor_locator1 = MaxNLocator(nbins=20)
    minor_locator2 = MaxNLocator(nbins=20)

    minor_locator1.set_params(integer=True)
    minor_locator2.set_params(integer=True)

    grid_helper = GridHelperCurveLinear(
        transform,
        extreme_finder=extreme_finder,
        grid_locator1=grid_locator1,
        grid_locator2=grid_locator2,
    )

    fig = plt.figure()

    ax = fig.add_subplot(1, 1, 1, axes_class=Axes, grid_helper=grid_helper)
    ax.set_aspect(aspect)
    ax.minorticks_on()

    trans = transform + ax.transData

    im = ax.pcolormesh(
        x,
        y,
        z,
        norm="linear",
        cmap="viridis",
        vmin=np.nanpercentile(z, pmin),
        vmax=np.nanpercentile(z, pmax),
        transform=trans,
        rasterized=True,
    )

    ax.set_xlabel(r"${}$".format(titles[0]))
    ax.set_ylabel(r"${}$".format(titles[1]))
    ax.set_title(r"${}=[{},{}]$".format(titles[2], *integration))

    cb = fig.colorbar(im, ax=ax)
    cb.ax.minorticks_on()
    cb.ax.set_ylabel(r"$I$ [arb. unit]")

    fig.show()


load_data(filename, ws)

integration = make_slice(ws, thickness, value, normal)

(x, y), z, e = extract_data(ws + "_slice")

T, aspect, titles = projection_transform(ws, y, normal)

make_plot(x, y, z, T, aspect, titles, integration, pmin, pmax)

Helper functions#

`load_data(filename, ws)`#

Load an MD workspace only if it does not already exist in the Mantid ADS.

param str filename:: Path to the input NeXus MD file.
param str ws:: Name of the output workspace in the ADS.

This avoids reloading large files when the workspace is already present.

`make_slice(ws, thickness, value, normal=[0, 0, 1])`#

Integrate the workspace over a finite thickness centered on a selected value along one reciprocal-space direction.

param str ws:: Name of the input workspace.
param float thickness:: Half-thickness of the integration region.
param float value:: Center of the slice.
param list normal:: Slice normal, expected to select one axis such as [1, 0, 0], [0, 1, 0], or [0, 0, 1].
return:: Lower and upper integration limits.
rtype:: list

The integration interval is constructed as [value - thickness, value + thickness]. The selected axis is passed to IntegrateMDHistoWorkspace through P1Bin, P2Bin, or P3Bin.

Note

This helper assumes the slice normal is aligned with one of the three MD dimensions. For arbitrary directions, a different approach would be needed.

`extract_data(ws)`#

Extract the non-integrated coordinates, signal, and errors from a sliced workspace.

param str ws:: Name of the sliced workspace.
return:: Meshgrid coordinates, signal array, and error array.
rtype:: tuple

The coordinate arrays are generated from the bin centers of the non-integrated workspace dimensions.

`projection_transform(ws, y, normal=[0, 0, 1])`#

Construct the affine transform used to display the reciprocal-space slice with projection-aware axes.

param str ws:: Name of the original workspace.
param numpy.ndarray y:: Second plotting coordinate array, used to anchor the affine offset.
param list normal:: Slice normal used to determine the integrated dimension.
return:: Affine transform matrix, aspect scaling, and axis titles.
rtype:: tuple

This helper:

retrieves the oriented lattice B matrix,
reads the W_MATRIX log from the run,
forms the projected basis B W,
uses QR and Cholesky factorizations to derive an in-plane transform, and
returns labels for the two plotted dimensions plus the integrated dimension.

`make_plot(x, y, z, T, aspect, titles, integration, pmin=5, pmax=95)`#

Plot the 2D slice using Matplotlib and axisartist.

param numpy.ndarray x:: First meshgrid coordinate array.
param numpy.ndarray y:: Second meshgrid coordinate array.
param numpy.ndarray z:: Signal values to plot.
param numpy.ndarray T:: 3 x 3 affine transform matrix.
param float aspect:: Aspect ratio returned by projection_transform.
param list titles:: Axis labels and integrated-dimension label.
param list integration:: Lower and upper integration limits.
param float pmin:: Lower percentile for color scaling.
param float pmax:: Upper percentile for color scaling.

The image is displayed with:

transformed coordinates using Affine2D,
a curved grid from GridHelperCurveLinear, and
percentile-based vmin and vmax to make contrast selection easier.

Typical usage#

For a slice normal to the third dimension:

load_data(filename, ws)
integration = make_slice(ws, thickness=0.1, value=0.0, normal=[0, 0, 1])
(x, y), z, e = extract_data(ws + "_slice")
T, aspect, titles = projection_transform(ws, y, normal=[0, 0, 1])
make_plot(x, y, z, T, aspect, titles, integration, pmin=0, pmax=95)

This produces a 2D reciprocal-space slice with transformed axes and a colorbar scaled between the selected signal percentiles.

Practical notes#

thickness is a half-width, not the full integrated width.
normal should contain exactly one entry equal to 1 and the others 0.
The workspace is sliced into a new workspace named ws + "_slice".
Percentile clipping is often more robust than fixed color limits when the data contains a few very intense pixels.
projection_transform expects the workspace to contain both an oriented lattice and a W_MATRIX log entry.

Plotting#

Overview#

Example script#

Helper functions#

load_data(filename, ws)#

make_slice(ws, thickness, value, normal=[0, 0, 1])#

extract_data(ws)#

projection_transform(ws, y, normal=[0, 0, 1])#

make_plot(x, y, z, T, aspect, titles, integration, pmin=5, pmax=95)#

Typical usage#

Practical notes#

`load_data(filename, ws)`#

`make_slice(ws, thickness, value, normal=[0, 0, 1])`#

`extract_data(ws)`#

`projection_transform(ws, y, normal=[0, 0, 1])`#

`make_plot(x, y, z, T, aspect, titles, integration, pmin=5, pmax=95)`#