Given our precomputed .datastack files present in our Nextcloud storage we can now visualize our MA-XRF data in a Dashboard viewer after downloading them once.
The total size of these 16 files is pretty big (56.5 GB) and downloading on my computer took about half an hour!
from fairdatanow import DataViewer
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import os
import numpy as np
import moseley as mos
import maxrf4u as mx
import re
import imageio as iio
from falnama import get_roi_imsconfiguration = {
'url': "https://laboppad.nl/falnama-project",
'user': os.getenv('NC_AUTH_USER'),
'password': os.getenv('NC_AUTH_PASS')
}filters = {'extensions': ['.datastack']}
dv = DataViewer(configuration, **filters)
datastack_files = dv.download_filtered()Ready with downloading 16 selected remote files to local cache: /home/frank/.cache/fairdatanow Ok, let’s start our exploration by plotting a thumbnail overview of all scanned areas.
imvis_list = []
extent_list = []
for datastack_file in datastack_files:
ds = mx.DataStack(datastack_file)
im = ds.read('imvis_reg')
extent = ds.read('imvis_extent')
imvis_list.append(im)
extent_list.append(extent)
object_nums = [re.sub(r'.*(WM-71803-\d\d).*', r'\1', dsf) for dsf in datastack_files]
fig, axs = plt.subplots(ncols=8, nrows=2, figsize=[15, 6])
axs = axs.flatten()
for i, [im, ax, extent] in enumerate(zip(imvis_list, axs, extent_list)):
ax.imshow(im, extent=extent)
ax.set_title(f'[{i}] {object_nums[i]}', fontsize=8)
ax.axis('off')
It came as a shock to me today that my two sulfur maps (i.e. a slice map and an NMF map) both are worse than the Bruker map.
filters = {'search': '.*Results/04.*71803-35_S.tif', 'use_regex': True}
dv = DataViewer(configuration, **filters)
bruker_sulfur_tif = dv.download_filtered()[0]Ready with downloading 1 selected remote files to local cache: /home/frank/.cache/fairdatanow bruker_S_map = plt.imread(bruker_sulfur_tif)fig, ax = plt.subplots()
ax.imshow(bruker_S_map)
ax.set_title('Bruker S map');
A first step now is to read relevant datasets from the datastack file and manually select regions of interest for page 35 (index 15)…
n = 15
datastack_file = datastack_files[n]
ds = mx.DataStack(datastack_file)
object_num = object_nums[n]
cube = ds.read('maxrf_cube', compute=False)
imvis_highres = ds.read('imvis_reg_highres')
imvis = ds.read('imvis_reg')
extent = ds.read('imvis_extent')
x_keVs = ds.read('maxrf_energies')
y_maxspectrum = ds.read('maxrf_maxspectrum')
nmf_elementmaps = ds.read('nmf_elementmaps')
#ppa = mx.Peak_Pattern_Atlas(datastack_file=datastack_file, tube_keV=23)
ppa = mos.PeakPatternAtlas(x_keVs=x_keVs)
ds.tree()Ready building Peak Pattern Atlas! / ├── compton_peak_energy (1,) float64 ├── hotmax_baselines (31, 4096) float64 ├── hotmax_noiselines (31, 4096) float64 ├── hotmax_peak_idxs_flat (34,) int64 ├── hotmax_peak_idxs_list (31, 2) int64 ├── hotmax_spectra (31, 4096) float32 ├── hotmax_spots (31, 2) int64 ├── hotmax_subpeak_idxs_list (31, 18) int64 ├── imvis_extent (4,) int64 ├── imvis_reg (617, 357, 3) uint8 ├── imvis_reg_highres (8256, 4777, 3) uint8 ├── maxrf_cube (617, 357, 4096) float32 ├── maxrf_energies (4096,) float64 ├── maxrf_maxspectrum (4096,) float32 ├── maxrf_sumspectrum (4096,) float64 ├── nmf_atomnums (18,) int64 ├── nmf_elementmaps (18, 617, 357) float32 ├── nmf_gausscomponents (35, 4096) float32 ├── nmf_peakmaps (35, 617, 357) float32 └── nmf_peaks2elements_matrix (18, 34) float32
/home/frank/.cache/fairdatanow/falnama-project/OPENDATA/maxrf/datastacks/WM-71803-35_400_500_50.datastack:
# broken...
# mx.get_instrument_pattern(datastack_file)fig, axs = plt.subplots(ncols=5, nrows=2, figsize=[20, 15]) axs = axs.flatten() for ax in axs: ax.imshow(imvis_highres, extent=extent)
xylims_list = [[int(ax.get_xlim()[0]), int(ax.get_xlim()[1]), int(ax.get_ylim()[0]), int(ax.get_ylim()[1])] for ax in axs]
xylims_list = [[29, 51, 106, 84], [39, 88, 136, 91], [154, 183, 242, 208], [275, 296, 249, 227], [40, 59, 289, 269],
[72, 124, 441, 397], [146, 168, 462, 443], [139, 161, 393, 370], [245, 268, 454, 433], [197, 219, 300, 278]]
roi_ims = get_roi_ims(datastack_file, xylims_list)fig, axs = plt.subplots(ncols=len(roi_ims)+1, figsize=[20, 4], squeeze=True)
edgecolors = ['violet', 'cyan', 'blue', 'red', 'white', 'orange', 'brown', 'green', 'maroon', 'black']
axs[0].imshow(imvis_highres, extent=extent)
for i, xylims in enumerate(xylims_list):
x0, x1, y0, y1 = xylims
axs[0].add_patch(Rectangle([x0, y1], x1 -x0, y0 - y1, linewidth=1, edgecolor=edgecolors[i], facecolor='none'))
for i, [ax, im] in enumerate(zip(axs[1:], roi_ims)):
ax.imshow(im)
ax.set_title(f'[{i}] {edgecolors[i]}')Here is the latest still very rough code for an ROI tabs dashboard.
plt.close('all')Ok, we now need some pretty advanced ROI hotmax functionality. Ehm, let’s build that here from scratch…
import falnama as fn###### BUILD DASHBOARD FIGURES CANVAS ################### r
roi_elements_list = [['Pb', 'K', 'Ca', 'Fe'], #0 We see clear peaks for K and Ca in the spectrum and slice maps why is Ca not showing clearly in the Ca (NMF) map?
['Pb_1', 'K', 'Ca', 'As_1', 'S'], #1
['Pb', 'Fe', 'Ca', 'Hg'], #2
['Pb', 'Pb_1', 'S'], #3
['Pb', 'Fe', 'Ca'], #4
['Pb', 'As', 'S', 'Fe'], #5 Here we see that the NMF maps are more clear the the alpha slices
['Au', 'Pb', 'As'], #6 Clearly gold on red lead. The Pb alpha slice map is confusing, no Hg present in the Au
['Pb', 'As', 'S', 'Au'], #7 The yellow robe is clearly orpiment. Here again the the overlapping Pb alpha with As is confusing
['As', 'S', 'Au', 'Pb'], #8 According to the NMF map there should be no lead present in the green robe
['Pb']] #9
roi_cubes = [cube[l:k, i:j].compute() for i, j, k, l in xylims_list]
roi_ims = get_roi_ims(datastack_file, xylims_list)roi_n = 3
ncols = len(roi_elements_list[roi_n])
hspace = 0.02
wspace = 0.02
xylims = xylims_list[roi_n]
l, r, b, t = xylims
w = r - l
h = b - t
roi_extent = [0, w, h, 0]
fig = plt.figure(layout='constrained', figsize=(15, 15))
topfigs, midfig, botfig = fig.subfigures(nrows=3, hspace=hspace, height_ratios=[2, 1, 2], squeeze=True)
topleftfig, topmidfig, toprightfig = topfigs.subfigures(ncols=3, wspace=wspace, width_ratios=[1, 1, 2], squeeze=True)
ax_00 = topleftfig.subplots()
ax_01 = topmidfig.subplots()
ax_020, ax_021 = toprightfig.subplots(nrows=2, squeeze=True, sharex=True)
spectra_axs = midfig.subplots(ncols=ncols, squeeze=True, sharex=True) #, sharey=True)
maps_axs = botfig.subplots(nrows=2, ncols=ncols, sharex=True, sharey=True)
########## TOP #################
# RGB overview and RGB ROI
ax_00.imshow(imvis_highres, extent=extent)
for i, xylims in enumerate(xylims_list):
x0, x1, y0, y1 = xylims
ax_00.add_patch(Rectangle([x0, y1], x1 -x0, y0 - y1, linewidth=1, edgecolor=edgecolors[i], facecolor='none'))
ax_00.annotate(f'ROI {i}', xy=[x0, y1], xytext=[0,5], textcoords='offset points', color='w', fontsize=8)
ax_00.set_title(object_num)
ax_01.imshow(roi_ims[roi_n], extent=roi_extent)
ax_01.set_title(f'ROI {roi_n}')
# ROI spectral plots and peek pattern atlas
# compute roi cube (core) maxspectrum
hotmax_idxs, spectra = get_hotmax_spectra(roi_cubes[roi_n])
roi_cube = roi_cubes[roi_n]
roi_h, roi_w, roi_d = roi_cube.shape
roi_flat = roi_cube.reshape([-1, roi_d])
roi_maxspectrum = roi_flat.max(axis=0) # actually this is the roi core max spectrum
roi_elements = roi_elements_list[roi_n]
ax_020.plot(x_keVs, roi_maxspectrum, label='ROI maxspectrum')
ax_020.plot(x_keVs, y_maxspectrum, color='r', alpha=0.2, label='full page maxspectrum')
ax_020.set_xlim([-1, 25])
ax_021.set_xlabel('energy [keV]')
#ax_020.legend() # need to fix blocking of view
# plot peak pattern atlas
select_elems = roi_elements_list[roi_n]
select_elems = [re.sub('_.*', '', elem) for elem in select_elems]
ppa.plot_atlas(ax=ax_021)
for i, element in enumerate(roi_elements):
# based on element now compute element slice map and related hotmax properties
slice_map, slice_idx = get_slice_map(datastack_file, element)
hx, hy, roi_x, roi_y, hotmax_spectrum = get_roi_slice_hotmax(datastack_file, slice_idx, xylims_list[roi_n])
# add hotmax locations to TOP image plots
ax_01.scatter(roi_x, roi_y)
# add slice line to spectral plots
#spectra_axs[i].axvline(x_keVs[slice_idx], color='r')
# add hotmax locations to slice maps
maps_axs[0, i].scatter(hx, hy, color='w', marker='+')
# add hotmax spectra and element shadows to mid figure
spectra_axs[i].plot(x_keVs, hotmax_spectrum)
spectra_axs[i].scatter(x_keVs[slice_idx], hotmax_spectrum[slice_idx], marker='s', edgecolor='r', facecolor='none')
#plot_element_shadow(element, x_keVs, flip=1, ax=spectra_axs[i], y=hotmax_spectrum)
spectra_axs[i].plot(x_keVs, y_maxspectrum, color='r', alpha=0.15)
spectra_axs[i].set_xlim([-1, 25])
spectra_axs[i].set_title(f'{element} (slice {x_keVs[slice_idx]:.1f} keV)', fontsize=8)
plot_element_lines(element, x_keVs, ax=spectra_axs[i], color='r', y=hotmax_spectrum)
# add NMF map for element
nmf_map = pick_nmf_elementmap(datastack_file, nmf_elementmaps, element)
maps_axs[0, i].imshow(slice_map)
maps_axs[0, i].set_title(f'{element} (slice {x_keVs[slice_idx]:.1f} keV)', fontsize=8)
maps_axs[1, i].imshow(nmf_map)
maps_axs[1, i].set_title(f'{element} (NMF)', fontsize=8)
# add ROI rectangle
x0, x1, y0, y1 = xylims_list[roi_n]
maps_axs[0, i].add_patch(Rectangle([x0, y1], x1 -x0, y0 - y1, linewidth=1, edgecolor=edgecolors[roi_n], facecolor='none'))
maps_axs[0, i].annotate(f'ROI {roi_n}', xy=[x0, y1], xytext=[0,5], textcoords='offset points', color='w', fontsize=10)
maps_axs[1, i].add_patch(Rectangle([x0, y1], x1 -x0, y0 - y1, linewidth=1, edgecolor=edgecolors[roi_n], facecolor='none'))
maps_axs[1, i].annotate(f'ROI {roi_n}', xy=[x0, y1], xytext=[0,5], textcoords='offset points', color='w', fontsize=10)Ok, Let’s extract a global hotmax spectrum for Pb_1
Pb_slice_map, Pb_slice_idx = get_slice_map(datastack_file, 'Pb_1')
yi, xi = np.argwhere(Pb_slice_map == Pb_slice_map.max())[0]
Pb_hotmax = cube[yi, xi, :].compute()fig, axs = plt.subplots(ncols=2, figsize=[15, 10], width_ratios=[1, 3])
ax0, ax1 = axs.flatten()
ax0.imshow(Pb_slice_map)
ax0.scatter(xi, yi, color='w', marker='+')
ax0.set_title('Pb_1')
ax1.fill_between(x_keVs, Pb_hotmax / Pb_hotmax.max(), color='r', label='Pb hotmax')
h_list = [0.0001, 0.001, 0.005, 0.01, 0.05, 0.1]
colors = cm.viridis(np.linspace(0, 1, len(h_list)))
for i, h_mm in enumerate(h_list):
Pb_xrf = mos.ElementXRF('Pb', h_mm=h_mm, excitation_energy_keV=30, x_keVs=x_keVs, std=0.02)
Pb_x, Pb_y = Pb_xrf.spectrum_xy
ax1.plot(Pb_x, 100*Pb_y, color=colors[i], label=f'Pb ({h_mm})')
ax1.legend()
ax1.set_xlim([0, 18]);I need to have access to the peak indexes…
Pb_xrf.get_pattern_dict(){'elem': 'Pb',
'atomic_number': 82,
'alpha_keV': 12.609790734147028,
'name': 'Lead',
'peaks_xy': array([[1.26097907e+01, 3.57595202e-02],
[1.05442700e+01, 3.31442037e-02],
[1.47655527e+01, 3.64028035e-03],
[1.51666247e+01, 1.29834767e-03],
[9.18062517e+00, 1.29593281e-03],
[2.37242804e+00, 5.87008945e-04],
[1.13464139e+01, 4.67320502e-04]]),
'alpha_escape_keV': 10.869790734147028}plot_roi_dashboard (datastack_file, xylims_list, roi_elements_list)
Create tabbed plot for roi spectra.
plot_roi_peak_patterns (roi_cube, prominence=2, ax=None)
get_roi_ims (datastack_file, xylims_list)
Get high resolution roi images for xylims_list.
get_slice_map (datastack_file, element)
*Compute simple single channel maps from datastack_file cube for element alpha channel.
Add ’_1’ to element string for beta peak slice map.
Returns: slice_map, slice_idx*
pick_nmf_elementmap (datastack_file, nmf_elementmaps, element)
Faster get nmf element map
get_nmf_map (datastack_file, element)
Read computed NMF element map for element from datastack.
plot_roi_maps (roi_cubes, xylims_list, roi_n, elements, y_maxspectrum=None, exc_keV=23)
*Plot roi alpha maps for list of elements.
If y_maxspectrum is specified normalize with maximum at channel index.*
get_roi_maps (roi_cube, elements, x_keVs, exc_keV=23)
*Slice roi cube at alpha channel index positions for each element in elements.
Returns: alpha_idxs, roi_maps, roi_vmax_list*
plot_element_lines (element, x_keVs, y=None, ax=None, exc_keV=23, color=None)
Add emission lines to spectral plot.
get_element_lines (element, x_keVs, exc_keV=23)
*Get sorted element peak idxs.
Returns: peak_idxs_sorted*
plot_hotmax_spectra (spectra, x_keVs, hotmax_idxs, elems=None, xlim=[-1, 22], elem_colors={'Pb': 'k', 'Au': 'orange', 'Fe': 'brown', 'Ca': 'g', 'K': 'b', 'As': 'magenta', 'Na': 'blue', 'Ni': 'red', 'Cu': 'green', 'Ti': 'blue', 'Zn': 'red'}, y_maxspectrum=None)
*Plot roi hotmax spectra and add emission lines for list of elements elems.
Returns: axs*
get_hotmax_spectra (roi_cube)
*Get hotmax spectra for roi_cube.
Returns: hotspot_idxs,spectra`*
plot_roi_peak_patterns (roi_cube, prominence=2, ax=None)
get_roi_slice_hotmax (datastack_file, slice_idx, xylims)
*Find full map and roi coordinates (x, y, roi_x, roi_y) of slice map maximum at channel slice_idx within ROI xylims.
Returns: x, y, roi_x, roi_y, hotmax_spectrum*