!pip install matplotlib s3fs "xarray[io]"
import warnings
import zarr
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.signal import stft
from zarr.errors import ZarrUserWarning
# The level-2 stores currently carry both Zarr v2 and v3 metadata, which zarr 3.x
# warns about on every read. Suppress that benign warning; it will disappear once
# the data is fully migrated to Zarr v3.
warnings.filterwarnings("ignore", category=ZarrUserWarning)
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[2], line 9
7 from matplotlib.colors import LogNorm
8 from scipy.signal import stft
----> 9 from zarr.errors import ZarrUserWarning
11 # The level-2 stores currently carry both Zarr v2 and v3 metadata, which zarr 3.x
12 # warns about on every read. Suppress that benign warning; it will disappear once
13 # the data is fully migrated to Zarr v3.
14 warnings.filterwarnings("ignore", category=ZarrUserWarning)
ImportError: cannot import name 'ZarrUserWarning' from 'zarr.errors' (/srv/fair-mast/.docs-venv/lib/python3.12/site-packages/zarr/errors.py)
Level 2 Data#
In this notebook we demonstrate the variety of different data available in the FAIR MAST dataset. In this example we are using level 2 MAST data, which includes cropping, interpolation, calibration, etc. of each signal, as well as mapping each diagnostic group try and follow IMAS naming convetions. The level 2 data is well-indexed data and follows the FAIR principles. Shots are also filtered using the plasma current to remove shots which were only used for testing, comissioning, machine calibration etc.
First we need to connect to the remote S3 storage bucket to access the data. Each shot from MAST is stored as a seperate Zarr file.
Using fsspec and xarray we can remotely read data directly over the web. In the example below we also turn on local file caching, allowing us to avoid reading over the network multiple times.
shot_id = 30421
endpoint_url = 'https://s3.echo.stfc.ac.uk'
url = f's3://mast/level2/shots/{shot_id}.zarr'
# Get a handle to the remote file
store = zarr.storage.FsspecStore.from_url(
url,
storage_options=dict(
protocol='simplecache',
target_protocol="s3",
cache_storage='.cache',
target_options=dict(
anon=True, endpoint_url=endpoint_url, asynchronous=True
)
)
)
Summary Profiles#
The summary group provides a collection of general physics quantities for an experiment.
profiles = xr.open_zarr(store, group='summary')
plot_1d_profiles(profiles)
profiles
Pulse Schedule#
profiles = xr.open_zarr(store, group='pulse_schedule')
fig, axes = plt.subplots(2, 1, figsize=(10, 5))
axes = axes.flatten()
profiles['i_plasma'].plot(x='time', ax=axes[0])
profiles['n_e_line'].plot(x='time', ax=axes[1])
for ax in axes:
ax.grid('on', alpha=0.5)
plt.tight_layout()
profiles
Magnetics#
Magnetic diagnostics for equilibrium identification and plasma shape control.
profiles = xr.open_zarr(store, group='magnetics')
fig, axes = plt.subplots(4, 3, figsize=(8, 10))
axes = axes.flatten()
profiles['b_field_pol_probe_ccbv_field'].plot.line(x='time', ax=axes[0], add_legend=False)
profiles['b_field_pol_probe_obv_field'].plot.line(x='time', ax=axes[1], add_legend=False)
profiles['b_field_pol_probe_obr_field'].plot.line(x='time', ax=axes[2], add_legend=False)
profiles['b_field_pol_probe_omv_voltage'].plot.line(x='time_mirnov', ax=axes[3], add_legend=False)
profiles['b_field_pol_probe_cc_field'].plot.line(x='time_mirnov', ax=axes[4], add_legend=False)
profiles['b_field_tor_probe_cc_field'].plot.line(x='time_mirnov', ax=axes[5], add_legend=False)
profiles['b_field_tor_probe_saddle_field'].plot.line(x='time_saddle', ax=axes[6], add_legend=False)
profiles['b_field_tor_probe_saddle_voltage'].plot.line(x='time_saddle', ax=axes[7], add_legend=False)
profiles['b_field_tor_probe_omaha_voltage'].plot.line(x='time_omaha', ax=axes[8], add_legend=False)
profiles['flux_loop_flux'].plot.line(x='time', ax=axes[9], add_legend=False)
profiles['ip'].plot.line(x='time', ax=axes[10], add_legend=False)
for ax in axes[:11]:
ax.grid('on', alpha=0.5)
# remove the unused 12th axis in the 4x3 grid
fig.delaxes(axes[11])
plt.tight_layout()
profiles
Looking at the spectrogram of one of the mirnov coils can show us information about the MHD modes. Here we see several mode instabilities occuring before the plasma is lost.
ds = profiles['b_field_pol_probe_omv_voltage'].isel(b_field_pol_probe_omv_channel=1)
# Parameters to limit the number of frequencies
nperseg = 2000 # Number of points per segment
nfft = 2000 # Number of FFT points
# Compute the Short-Time Fourier Transform (STFT)
sample_rate = 1/(ds.time_mirnov[1] - ds.time_mirnov[0])
f, t, Zxx = stft(ds, fs=int(sample_rate), nperseg=nperseg, nfft=nfft)
fig, ax = plt.subplots(figsize=(15, 5))
cax = ax.pcolormesh(t, f/1000, np.abs(Zxx), shading='nearest', cmap='jet', norm=LogNorm(vmin=1e-5))
ax.set_ylim(0, 50)
ax.set_title(f'Shot {shot_id}')
ax.set_ylabel('Frequency [Hz]')
ax.set_xlabel('Time [sec]')
plt.colorbar(cax, ax=ax)
plt.show()
Spectrometer Visible#
Spectrometer in visible light range diagnostic
profiles = xr.open_zarr(store, group='spectrometer_visible')
profiles['filter_spectrometer_dalpha_voltage'].plot.line(x='time')
profiles['filter_spectrometer_bes_voltage'].isel(bes_channel=0).plot.line(x='time_bes')
profiles
PF Active#
profiles = xr.open_zarr(store, group='pf_active')
fig, axes = plt.subplots(2, 1, figsize=(10, 6))
# coil_current includes the central solenoid (the 'SOL' channel), which
# carries the largest current; the legend identifies each coil channel.
profiles['coil_current'].plot.line(x='time', ax=axes[0])
profiles['coil_voltage'].plot.line(x='time', ax=axes[1])
plt.tight_layout()
profiles
Soft X-rays#
profiles = xr.open_zarr(store, group='soft_x_rays')
fig, axes = plt.subplots(3, 1)
profiles['horizontal_cam_lower'].plot.line(x='time', ax=axes[1], add_legend=False)
axes[1].set_ylim(0, 0.02)
profiles['horizontal_cam_upper'].plot.line(x='time', ax=axes[2], add_legend=False)
axes[2].set_ylim(0, 0.02)
if "tangential_cam" in profiles:
profiles['tangential_cam'].plot.line(x='time', ax=axes[0], add_legend=False)
axes[0].set_ylim(0, 0.2)
plt.tight_layout()
profiles
Thomson Profiles#
Thomson scattering measurements in a tokamak provide information about the plasma’s electron temperature and density profiles. The diagnostic analyses the scattering of laser light off free electrons in the plasma from a number of radial channels.
Below we plot the following profiles measured by the Thomson diagnostic
\(T_e\) - Electron temperature
\(N_e\) - Electron density
\(P_e\) - Electron pressure
profiles = xr.open_zarr(store, group='thomson_scattering')
profiles
fig, axes = plt.subplots(3, 1)
axes = axes.flatten()
profiles.t_e.plot(x='time', y='major_radius', ax=axes[0])
profiles.n_e.plot(x='time', y='major_radius', ax=axes[1])
profiles.p_e.plot(x='time', y='major_radius', ax=axes[2])
plt.tight_layout()
profiles
fig, axes = plt.subplots(2, 1)
profiles['t_e_core'].plot(x='time', ax=axes[0])
profiles['n_e_core'].plot(x='time', ax=axes[1])
for ax in axes:
ax.grid('on', alpha=0.5)
plt.tight_layout()
profiles
CXRS Profiles#
Charge Exchange Recombination Spectroscopy (CXRS) measurements provide information about ion temperature and plasma rotation. This diagnostic analyses the light emitted from charge exchange reactions between injected neutral beams and plasma ions.
Below we plot the following profiles measured by the CXRS diagnostic
\(T_i\) - Ion temperature
\(V_i\) - Ion velocity
profiles = xr.open_zarr(store, group='charge_exchange')
fig, axes = plt.subplots(2, 1)
for name, ax in [('t_i', axes[0]), ('v_i', axes[1])]:
da = profiles[name]
da.plot(x='time', y='major_radius', ax=ax, vmax=1000)
# clip the axes to the region where the diagnostic actually has data
valid = da.notnull()
times = da['time'].where(valid.any('major_radius'), drop=True)
radii = da['major_radius'].where(valid.any('time'), drop=True)
ax.set_xlim(float(times.min()), float(times.max()))
ax.set_ylim(float(radii.min()), float(radii.max()))
plt.tight_layout()
profiles
Equilibrium#
profiles = xr.open_zarr(store, group='equilibrium')
profile_1d = profiles[["beta_tor_normal", "wmhd", "li", "elongation", "triangularity_upper", "q95", "vloop_dynamic", "ip_rating", "lcfs_r", "lcfs_z"]]
plot_1d_profiles(profile_1d)
profiles
fig, axes = plt.subplots(1, 3, figsize=(10, 5))
profiles['j_phi'].isel(time=50).plot(ax=axes[0], x='major_radius')
profiles['psi'].isel(time=50).plot(ax=axes[1], x='major_radius')
profiles['q'].isel(time=50).plot(ax=axes[2])
plt.tight_layout()
Gas Injection#
profiles = xr.open_zarr(store, group='gas_injection')
plot_1d_profiles(profiles[["inboard_total", "outboard_total", "pressure", "total_injected"]])
plt.tight_layout()
profiles