Introduction to Scipp

Multi-dimensional arrays with labeled dimensions and physical units

scipp.github.io

Scipp is an open-source library developed by ESS for handling, manipulating and visualizing multi-dimensional data arrays.

It enriches raw NumPy-like arrays by adding named dimensions and associated coordinates. In addition, it supports

• Physical units which are handled in arithmetic operations

• Histograms, i.e., bin-edge axes, which are by 1 longer than the data extent

• Propagation of uncertainties

%matplotlib inline
import numpy as np
import scipp as sc
import plopp as pp
import matplotlib.pyplot as plt

from scipp_utils import quiz, plot, scatter, fetch_data

rng = np.random.default_rng(seed=1234)

1. Labeled dimensions: why do we need them?

Say we have a 2D rectangular array of data

ny, nx = 10, 20
a = np.sin(np.arange(ny) / (ny / 4)).reshape((-1, 1)) * np.cos(np.arange(nx) / (ny / 4))
a.shape

(10, 20)

that looks like

plot(a)

The task is now to slice out row number 4. Because of the shape of the array, we know that the row dimension is the smallest, so we slice the first dimension of the 2D array:

# Slice out row number 4
plot(a[4, :])

We can’t always deduce from the shape

Now say we have an array which has a square shape:

ny, nx = 20, 20
a = np.sin(np.arange(ny) / (ny / 4)).reshape((-1, 1)) * np.cos(np.arange(nx) / (ny / 4))
a.shape

(20, 20)

plot(a)

Do we slice the first or the second index of the 2D array?

# Not always obvious which dimension is which
plot(a[:, 4], a[4, :])

The situation gets worse with more dimensions

Say we now have an array that has 4 dimensions: x, y, z, t (in that order, maybe?, or is it z, y, x, t, or t, x, y, z?)

a = np.random.random([20] * 4)
a.shape

(20, 20, 20, 20)

Quiz time!

quiz(1)

Introducing labeled dimensions

        

Xarray introduced labels to multi-dimensional Numpy arrays.

“real-world datasets are usually more than just raw numbers; they have labels which encode information about how the array values map to locations in space, time, etc.”

We have embraced, and to a large extent copied, the Xarray mechanism.

var = sc.array(dims=["x", "y", "z", "time"], values=a)
var

scipp.Variable (1.22 MB)

• • (x: 20, y: 20, z: 20, time: 20)
    float64
    𝟙
    0.927, 0.500, ..., 0.728, 0.234

    Values:
    array([[[[0.92724983, 0.49993413, 0.16599067, ..., 0.91301157,
              0.35727723, 0.98606653],
             [0.89467734, 0.36372609, 0.94366775, ..., 0.64975446,
              0.8487413 , 0.83129484],
             [0.20655218, 0.78091904, 0.06175176, ..., 0.39330951,
              0.99274471, 0.32103939],
             ...,
             [0.68342668, 0.84232305, 0.53420582, ..., 0.365706  ,
              0.76378493, 0.18323771],
             [0.7767901 , 0.19251071, 0.06738654, ..., 0.69425547,
              0.24296738, 0.14260514],
             [0.43267841, 0.89532852, 0.25946203, ..., 0.64006416,
              0.45186943, 0.65812809]],
    
            [[0.94393598, 0.22677644, 0.6942398 , ..., 0.60914589,
              0.7455753 , 0.45471123],
             [0.06697583, 0.50750317, 0.04548286, ..., 0.80627193,
              0.67058966, 0.45525376],
             [0.63391356, 0.43284414, 0.83211742, ..., 0.27118052,
              0.71375232, 0.61020412],
             ...,
             [0.73091467, 0.71810253, 0.84506733, ..., 0.74241337,
              0.27912304, 0.6055827 ],
             [0.49393801, 0.0923407 , 0.13064871, ..., 0.44513409,
              0.94900911, 0.66972311],
             [0.24875631, 0.82255538, 0.34421695, ..., 0.61092445,
              0.24005149, 0.94670113]],
    
            [[0.89832867, 0.43421003, 0.45586941, ..., 0.2253544 ,
              0.71787842, 0.76321177],
             [0.49232674, 0.08450241, 0.76131766, ..., 0.55686949,
              0.76884559, 0.73979538],
             [0.83860711, 0.02719219, 0.5819619 , ..., 0.62909381,
              0.64307022, 0.54670259],
             ...,
             [0.46096215, 0.55724201, 0.59020593, ..., 0.93257896,
              0.95277287, 0.36232569],
             [0.96773279, 0.00175159, 0.75546769, ..., 0.09265446,
              0.2697582 , 0.21847838],
             [0.97484402, 0.98512866, 0.83415194, ..., 0.36129847,
              0.01342972, 0.13023939]],
    
            ...,
    
            [[0.72173368, 0.51042845, 0.51698085, ..., 0.53353585,
              0.73913658, 0.37643072],
             [0.87439392, 0.40945752, 0.36641432, ..., 0.53628975,
              0.56901655, 0.81635216],
             [0.10972303, 0.43785125, 0.75249457, ..., 0.91799771,
              0.29046061, 0.85723201],
             ...,
             [0.48491589, 0.66231988, 0.00905441, ..., 0.14521791,
              0.93529346, 0.31020824],
             [0.32488816, 0.52412851, 0.18011863, ..., 0.09867884,
              0.93879242, 0.49164102],
             [0.94796108, 0.1750707 , 0.22577889, ..., 0.12074863,
              0.53968286, 0.57447579]],
    
            [[0.76508451, 0.36291907, 0.58293725, ..., 0.6251778 ,
              0.87575522, 0.36321896],
             [0.9274272 , 0.18045275, 0.50503968, ..., 0.39246741,
              0.91011043, 0.38328661],
             [0.88931734, 0.65049502, 0.40432287, ..., 0.25047324,
              0.43599181, 0.71133785],
             ...,
             [0.16308953, 0.36933598, 0.96580477, ..., 0.56171509,
              0.72598289, 0.16826535],
             [0.75675638, 0.75247268, 0.81410246, ..., 0.07352113,
              0.35109199, 0.92102123],
             [0.7295137 , 0.26286687, 0.7556685 , ..., 0.30413782,
              0.54940719, 0.29466013]],
    
            [[0.27454191, 0.70967315, 0.18512923, ..., 0.94616267,
              0.05073672, 0.36302733],
             [0.64579798, 0.88497406, 0.92423522, ..., 0.97919066,
              0.45926433, 0.23708913],
             [0.81979253, 0.45358186, 0.16817827, ..., 0.03268104,
              0.35256736, 0.7862389 ],
             ...,
             [0.26297646, 0.67212519, 0.70910651, ..., 0.11937525,
              0.79847872, 0.27576807],
             [0.71845744, 0.11836819, 0.29526798, ..., 0.25825722,
              0.04899791, 0.81934253],
             [0.58522082, 0.89718135, 0.36649658, ..., 0.61517446,
              0.99636598, 0.62297896]]],
    
    
           [[[0.08624627, 0.18850385, 0.19003988, ..., 0.54660584,
              0.85278244, 0.79396748],
             [0.28982592, 0.68410386, 0.85358409, ..., 0.48313971,
              0.83366813, 0.0336281 ],
             [0.15687452, 0.05723389, 0.88132811, ..., 0.38271037,
              0.87274314, 0.89375888],
             ...,
             [0.2537066 , 0.09473751, 0.6548335 , ..., 0.05740774,
              0.88719714, 0.62679162],
             [0.96112388, 0.21231   , 0.13628953, ..., 0.89287046,
              0.8792287 , 0.64389663],
             [0.79483039, 0.38352238, 0.55240888, ..., 0.10581238,
              0.01814151, 0.1632724 ]],
    
            [[0.66281628, 0.95858813, 0.93156885, ..., 0.99381714,
              0.92303506, 0.48901655],
             [0.92508841, 0.19828963, 0.90706194, ..., 0.85538948,
              0.06484454, 0.61265777],
             [0.93634635, 0.87976573, 0.57118959, ..., 0.452768  ,
              0.32563774, 0.69984692],
             ...,
             [0.23210335, 0.94652703, 0.01097254, ..., 0.65636539,
              0.7545695 , 0.56657687],
             [0.445445  , 0.19782987, 0.54197171, ..., 0.32304741,
              0.56668897, 0.60236579],
             [0.11993505, 0.14655708, 0.77407612, ..., 0.06179385,
              0.94868283, 0.97470606]],
    
            [[0.25598675, 0.46251897, 0.8818316 , ..., 0.34614642,
              0.23471103, 0.74820022],
             [0.07779633, 0.00921728, 0.47026266, ..., 0.27199877,
              0.83855086, 0.65529   ],
             [0.68672776, 0.37478465, 0.89913834, ..., 0.90940165,
              0.87749016, 0.9862704 ],
             ...,
             [0.11348508, 0.03565039, 0.61534012, ..., 0.72248969,
              0.89444634, 0.92029683],
             [0.36931007, 0.18851492, 0.3636989 , ..., 0.70154413,
              0.0589707 , 0.26837041],
             [0.14167886, 0.37724852, 0.22616596, ..., 0.16727586,
              0.33061671, 0.28517491]],
    
            ...,
    
            [[0.62417366, 0.81876866, 0.85936019, ..., 0.91697016,
              0.36785912, 0.6800574 ],
             [0.88203089, 0.64089518, 0.94314552, ..., 0.39490175,
              0.75122327, 0.53543667],
             [0.47366449, 0.38843916, 0.22863509, ..., 0.67205907,
              0.5087258 , 0.71551018],
             ...,
             [0.15367012, 0.32893913, 0.21670339, ..., 0.44182116,
              0.77695322, 0.746167  ],
             [0.65665677, 0.4354982 , 0.8416897 , ..., 0.51926931,
              0.92298771, 0.39851988],
             [0.79501885, 0.67148711, 0.87514043, ..., 0.22341338,
              0.51473346, 0.06678256]],
    
            [[0.45046213, 0.10274093, 0.53244226, ..., 0.23473611,
              0.28015792, 0.64784495],
             [0.15659579, 0.41962717, 0.34700396, ..., 0.87198257,
              0.64633263, 0.87547797],
             [0.44938402, 0.01764845, 0.90154034, ..., 0.57836837,
              0.91548124, 0.05431431],
             ...,
             [0.17312099, 0.76457741, 0.50035832, ..., 0.40349914,
              0.65687191, 0.95519872],
             [0.36195193, 0.73767074, 0.73801568, ..., 0.29758198,
              0.42054157, 0.86045471],
             [0.81373922, 0.68995423, 0.98693446, ..., 0.63083565,
              0.68361094, 0.602181  ]],
    
            [[0.97463268, 0.95949589, 0.21287622, ..., 0.74782125,
              0.64890746, 0.56479703],
             [0.02156447, 0.03873882, 0.56308736, ..., 0.28958565,
              0.03353361, 0.6891188 ],
             [0.73884044, 0.47666106, 0.95293318, ..., 0.28185597,
              0.79224453, 0.56377956],
             ...,
             [0.55315023, 0.32554078, 0.74159489, ..., 0.44018017,
              0.01768048, 0.52417292],
             [0.69907602, 0.15303219, 0.19535103, ..., 0.59091888,
              0.30785111, 0.53494023],
             [0.07676794, 0.4734298 , 0.9612888 , ..., 0.39236042,
              0.19205521, 0.57726093]]],
    
    
           [[[0.00479158, 0.68425181, 0.88348308, ..., 0.30122089,
              0.29450532, 0.81268838],
             [0.19041542, 0.08886042, 0.99846876, ..., 0.41917396,
              0.00706894, 0.60245423],
             [0.81805515, 0.39875972, 0.85211298, ..., 0.81366636,
              0.4417271 , 0.61593453],
             ...,
             [0.85226375, 0.78200751, 0.07214408, ..., 0.45229157,
              0.98711909, 0.01820432],
             [0.71450869, 0.45884559, 0.71697193, ..., 0.41742177,
              0.15044197, 0.58229837],
             [0.84885655, 0.66564493, 0.46663477, ..., 0.0047362 ,
              0.30462879, 0.21739544]],
    
            [[0.98869749, 0.41658116, 0.61677267, ..., 0.84888609,
              0.73774736, 0.71167644],
             [0.95844391, 0.20524544, 0.38600718, ..., 0.6113765 ,
              0.76219409, 0.31119932],
             [0.64239919, 0.73728704, 0.53909083, ..., 0.54912691,
              0.59228143, 0.2750154 ],
             ...,
             [0.73352224, 0.29462478, 0.25150342, ..., 0.96595329,
              0.98972797, 0.7726631 ],
             [0.30404941, 0.68891998, 0.49549254, ..., 0.09536449,
              0.54956649, 0.01256276],
             [0.3675224 , 0.37741549, 0.28184851, ..., 0.83220656,
              0.3361944 , 0.0337611 ]],
    
            [[0.4259418 , 0.57015747, 0.76956419, ..., 0.25738308,
              0.24177076, 0.88946012],
             [0.67726954, 0.32597229, 0.10201559, ..., 0.0318805 ,
              0.21083492, 0.59638205],
             [0.31132538, 0.11060485, 0.57001942, ..., 0.07073551,
              0.41619257, 0.34670155],
             ...,
             [0.97521759, 0.8138224 , 0.47792506, ..., 0.68429268,
              0.92418903, 0.42147159],
             [0.14115954, 0.74216744, 0.37217945, ..., 0.31566865,
              0.54838792, 0.07084265],
             [0.4217462 , 0.31732775, 0.10334198, ..., 0.60973643,
              0.0844357 , 0.77785849]],
    
            ...,
    
            [[0.72645395, 0.01869636, 0.34012839, ..., 0.16751079,
              0.35484258, 0.93526113],
             [0.99212867, 0.79079135, 0.57811128, ..., 0.97369116,
              0.71305513, 0.62819495],
             [0.67901957, 0.01362673, 0.32241452, ..., 0.6551284 ,
              0.40713906, 0.13126122],
             ...,
             [0.89295328, 0.14609531, 0.68799659, ..., 0.09961977,
              0.06422657, 0.16153391],
             [0.73633967, 0.72602473, 0.67980452, ..., 0.72135871,
              0.25105185, 0.90707738],
             [0.28572609, 0.96807734, 0.9521234 , ..., 0.38788493,
              0.59288325, 0.26974114]],
    
            [[0.7422356 , 0.69821868, 0.01050593, ..., 0.83100302,
              0.69784297, 0.11883557],
             [0.10743448, 0.95528608, 0.10049781, ..., 0.62775872,
              0.44861185, 0.98367752],
             [0.52600783, 0.35873139, 0.62186456, ..., 0.73044626,
              0.57083943, 0.73342398],
             ...,
             [0.10048095, 0.72439759, 0.58455655, ..., 0.88404982,
              0.34778152, 0.24023127],
             [0.0740435 , 0.87093244, 0.27403956, ..., 0.71591895,
              0.54344372, 0.00735904],
             [0.80827644, 0.26241667, 0.56521661, ..., 0.17469006,
              0.95559996, 0.25296523]],
    
            [[0.05066376, 0.00233922, 0.91602745, ..., 0.54545342,
              0.38248383, 0.65058518],
             [0.8798513 , 0.42807066, 0.20463774, ..., 0.98640514,
              0.65110079, 0.15640891],
             [0.85671568, 0.94645848, 0.52901419, ..., 0.1543499 ,
              0.15412193, 0.06081256],
             ...,
             [0.49306436, 0.49661425, 0.92486937, ..., 0.43218326,
              0.3384825 , 0.69921014],
             [0.47978898, 0.60898682, 0.30445754, ..., 0.30174054,
              0.94060035, 0.82121712],
             [0.73455118, 0.95812571, 0.18745742, ..., 0.48201067,
              0.37181443, 0.80422533]]],
    
    
           ...,
    
    
           [[[0.75444435, 0.21052252, 0.25473514, ..., 0.56793984,
              0.84048457, 0.07904822],
             [0.54529911, 0.32666102, 0.06591715, ..., 0.88894725,
              0.12185084, 0.18035185],
             [0.27416782, 0.10061743, 0.93586322, ..., 0.75854822,
              0.90477422, 0.40402628],
             ...,
             [0.52440837, 0.44059528, 0.27088888, ..., 0.56036428,
              0.69749182, 0.3205202 ],
             [0.88851907, 0.02978186, 0.28102694, ..., 0.68030818,
              0.97268153, 0.81931492],
             [0.94878257, 0.44999111, 0.22321753, ..., 0.64618304,
              0.7940522 , 0.15992511]],
    
            [[0.68494519, 0.43874286, 0.43650904, ..., 0.63124074,
              0.40689382, 0.53424961],
             [0.39801839, 0.20483552, 0.02522938, ..., 0.65606512,
              0.54773804, 0.10911548],
             [0.62991986, 0.02212523, 0.56262157, ..., 0.44648921,
              0.94205421, 0.8013069 ],
             ...,
             [0.46941601, 0.77045896, 0.62581482, ..., 0.83156907,
              0.27774359, 0.22916431],
             [0.11698349, 0.82148261, 0.52888626, ..., 0.88924832,
              0.02084974, 0.59472427],
             [0.54930151, 0.91420514, 0.43518764, ..., 0.28039537,
              0.7937925 , 0.70757599]],
    
            [[0.01409606, 0.37321952, 0.37048782, ..., 0.42755429,
              0.71710973, 0.13540134],
             [0.76252487, 0.46334547, 0.38577945, ..., 0.93194833,
              0.6794517 , 0.06092657],
             [0.86090196, 0.52616054, 0.612498  , ..., 0.48965827,
              0.19763517, 0.67498411],
             ...,
             [0.38054892, 0.62800535, 0.32890727, ..., 0.54331461,
              0.07324477, 0.51980274],
             [0.5796129 , 0.1440572 , 0.34206927, ..., 0.15563632,
              0.16856375, 0.66705702],
             [0.14403946, 0.70696419, 0.476766  , ..., 0.74935442,
              0.58517535, 0.73553549]],
    
            ...,
    
            [[0.07041195, 0.67529144, 0.13143863, ..., 0.10048729,
              0.13277627, 0.935807  ],
             [0.45724803, 0.62867682, 0.04905815, ..., 0.26508899,
              0.59492088, 0.88039491],
             [0.03391257, 0.36092076, 0.84801561, ..., 0.83260567,
              0.83638589, 0.61859684],
             ...,
             [0.17351119, 0.40988219, 0.98896572, ..., 0.2272736 ,
              0.52301804, 0.76062681],
             [0.35204221, 0.98983358, 0.15583589, ..., 0.27401562,
              0.52182627, 0.65780499],
             [0.86478474, 0.83605172, 0.89176497, ..., 0.35258153,
              0.07196896, 0.96125879]],
    
            [[0.13980907, 0.7875445 , 0.61815414, ..., 0.1451625 ,
              0.8627399 , 0.13818134],
             [0.14002531, 0.13281348, 0.00172431, ..., 0.47776522,
              0.26794415, 0.19479731],
             [0.0518809 , 0.13966449, 0.0621094 , ..., 0.67657577,
              0.06950385, 0.6456786 ],
             ...,
             [0.77379764, 0.48532798, 0.63447915, ..., 0.32032413,
              0.81410313, 0.55822408],
             [0.54176783, 0.16325825, 0.84903648, ..., 0.76579073,
              0.6833354 , 0.59774696],
             [0.83268427, 0.8755876 , 0.9297804 , ..., 0.65751921,
              0.77548648, 0.5574279 ]],
    
            [[0.46865885, 0.27220997, 0.3580191 , ..., 0.35170523,
              0.58851203, 0.91794438],
             [0.92067471, 0.37195063, 0.14243789, ..., 0.48244798,
              0.04263322, 0.41736029],
             [0.55241872, 0.99267812, 0.73130585, ..., 0.25177584,
              0.95652357, 0.7645796 ],
             ...,
             [0.41108404, 0.86019141, 0.57644868, ..., 0.59295109,
              0.46224371, 0.11727564],
             [0.53517694, 0.16154263, 0.41073609, ..., 0.45248425,
              0.01091592, 0.00152415],
             [0.86136255, 0.70413993, 0.36409386, ..., 0.1102886 ,
              0.18630199, 0.50953568]]],
    
    
           [[[0.85712706, 0.00593189, 0.23561811, ..., 0.72605492,
              0.42190059, 0.5569326 ],
             [0.21694134, 0.32093344, 0.60189147, ..., 0.93541231,
              0.98807067, 0.87076622],
             [0.36589059, 0.99839984, 0.26525606, ..., 0.05262929,
              0.11604961, 0.99806939],
             ...,
             [0.87152851, 0.16598927, 0.53272232, ..., 0.25311283,
              0.62511547, 0.80986077],
             [0.19950667, 0.31243649, 0.92409257, ..., 0.50597519,
              0.19424482, 0.87929482],
             [0.37355987, 0.77574002, 0.80522809, ..., 0.16696746,
              0.07111236, 0.7514431 ]],
    
            [[0.76063032, 0.25689224, 0.37195464, ..., 0.54598555,
              0.8924825 , 0.63199575],
             [0.59216855, 0.89553363, 0.98876219, ..., 0.62691885,
              0.49396663, 0.67699688],
             [0.82949371, 0.80227841, 0.0725771 , ..., 0.85926359,
              0.67092379, 0.3868004 ],
             ...,
             [0.94122718, 0.821609  , 0.14575985, ..., 0.14521755,
              0.28683765, 0.90894281],
             [0.07087918, 0.54518341, 0.81640691, ..., 0.17952289,
              0.17804276, 0.32758871],
             [0.05809748, 0.02098363, 0.83001556, ..., 0.33630461,
              0.38202689, 0.19450411]],
    
            [[0.08202932, 0.6902656 , 0.17096062, ..., 0.76829636,
              0.02772487, 0.17865257],
             [0.73186754, 0.94768927, 0.15566548, ..., 0.15165346,
              0.20575034, 0.1207281 ],
             [0.92086162, 0.0666654 , 0.12493554, ..., 0.9615467 ,
              0.63644258, 0.58632348],
             ...,
             [0.47751801, 0.39630213, 0.44738085, ..., 0.42492748,
              0.25802237, 0.8247561 ],
             [0.26242569, 0.29247922, 0.77817143, ..., 0.52193554,
              0.48870241, 0.84501297],
             [0.39192758, 0.02957993, 0.27292987, ..., 0.96842536,
              0.00576467, 0.21963449]],
    
            ...,
    
            [[0.45400832, 0.81076218, 0.21935231, ..., 0.43303348,
              0.2692096 , 0.00846294],
             [0.84716239, 0.12961095, 0.71205685, ..., 0.67280949,
              0.81242063, 0.75608572],
             [0.63915091, 0.56080803, 0.80969547, ..., 0.54103455,
              0.56482727, 0.03392996],
             ...,
             [0.80082975, 0.02783454, 0.591219  , ..., 0.27252893,
              0.39541249, 0.91446232],
             [0.24857981, 0.65095977, 0.28553415, ..., 0.71660406,
              0.40218732, 0.61835936],
             [0.16255042, 0.33846003, 0.47295493, ..., 0.20514029,
              0.72638479, 0.05847047]],
    
            [[0.53541442, 0.42371501, 0.74779775, ..., 0.67637444,
              0.3463082 , 0.7729033 ],
             [0.04146092, 0.14861115, 0.31186091, ..., 0.54460686,
              0.30590373, 0.50805071],
             [0.24912934, 0.55783867, 0.28827699, ..., 0.35315835,
              0.53196954, 0.97545888],
             ...,
             [0.69018968, 0.21620904, 0.90802715, ..., 0.24737809,
              0.41714848, 0.87164634],
             [0.23568319, 0.25884093, 0.64340974, ..., 0.65834343,
              0.12914308, 0.06366576],
             [0.51368173, 0.93095775, 0.63710151, ..., 0.40176416,
              0.79063637, 0.10059557]],
    
            [[0.69651702, 0.546848  , 0.71198762, ..., 0.86963897,
              0.6694778 , 0.35403408],
             [0.1576718 , 0.37828942, 0.37826664, ..., 0.08005755,
              0.40252224, 0.694326  ],
             [0.78455499, 0.2053314 , 0.52520319, ..., 0.13049379,
              0.6555466 , 0.56545245],
             ...,
             [0.99104932, 0.76010849, 0.3039156 , ..., 0.95195146,
              0.81122065, 0.63320369],
             [0.6510016 , 0.16638467, 0.68727973, ..., 0.92788047,
              0.51861186, 0.1102183 ],
             [0.10044871, 0.02422699, 0.17981271, ..., 0.72447215,
              0.08893293, 0.54374532]]],
    
    
           [[[0.54876883, 0.97517435, 0.53393699, ..., 0.75528858,
              0.86123632, 0.71484477],
             [0.98795194, 0.92821684, 0.57089299, ..., 0.99166329,
              0.07981442, 0.0779175 ],
             [0.12448897, 0.87246789, 0.67619932, ..., 0.93270259,
              0.5957397 , 0.86602232],
             ...,
             [0.22620057, 0.14288701, 0.53249006, ..., 0.05197059,
              0.77473409, 0.57808495],
             [0.94427023, 0.86223033, 0.31368485, ..., 0.32841413,
              0.32851205, 0.4878051 ],
             [0.36893959, 0.17589367, 0.99961457, ..., 0.96164051,
              0.36553845, 0.16459906]],
    
            [[0.30821069, 0.01032365, 0.09374503, ..., 0.9294343 ,
              0.54407655, 0.61835432],
             [0.71275477, 0.2198664 , 0.14529723, ..., 0.41192777,
              0.2507902 , 0.43312639],
             [0.59520607, 0.38989215, 0.10056406, ..., 0.37694754,
              0.37175386, 0.85390237],
             ...,
             [0.6012846 , 0.11002837, 0.7555092 , ..., 0.2313444 ,
              0.11670829, 0.24144134],
             [0.5679061 , 0.61446513, 0.21898892, ..., 0.40847922,
              0.18044447, 0.82917792],
             [0.87000224, 0.45107913, 0.24591306, ..., 0.69777113,
              0.64221148, 0.1104425 ]],
    
            [[0.58644036, 0.97173298, 0.15911906, ..., 0.37790684,
              0.88080998, 0.95715579],
             [0.74600631, 0.36200674, 0.63350868, ..., 0.64614041,
              0.73664039, 0.29961865],
             [0.15225974, 0.34229525, 0.29609715, ..., 0.77518648,
              0.20050935, 0.85067728],
             ...,
             [0.96976075, 0.29038724, 0.84606013, ..., 0.29156929,
              0.59113964, 0.03738946],
             [0.90753178, 0.69903723, 0.18622374, ..., 0.59349386,
              0.84321672, 0.04028122],
             [0.08630001, 0.70963556, 0.4768245 , ..., 0.475334  ,
              0.367331  , 0.32079771]],
    
            ...,
    
            [[0.20470871, 0.61996625, 0.07614679, ..., 0.72119104,
              0.64088126, 0.17990009],
             [0.07492384, 0.46218628, 0.95323946, ..., 0.26682939,
              0.8719155 , 0.08553354],
             [0.31571407, 0.17070284, 0.0878798 , ..., 0.95404076,
              0.13351047, 0.82134284],
             ...,
             [0.5738649 , 0.45368918, 0.6700238 , ..., 0.36009561,
              0.12829481, 0.34562304],
             [0.59252535, 0.28065799, 0.58755818, ..., 0.09276664,
              0.11757979, 0.26335586],
             [0.01664826, 0.95859783, 0.61540261, ..., 0.50998329,
              0.47590321, 0.43848769]],
    
            [[0.42024112, 0.78320008, 0.37225494, ..., 0.31050587,
              0.65967409, 0.99779411],
             [0.44347361, 0.09245578, 0.68340645, ..., 0.07214126,
              0.72799055, 0.53754808],
             [0.90815038, 0.20598156, 0.45891336, ..., 0.816346  ,
              0.00314905, 0.82672805],
             ...,
             [0.437584  , 0.45642369, 0.41356332, ..., 0.71533411,
              0.82451311, 0.25067518],
             [0.0924358 , 0.10875517, 0.1414318 , ..., 0.75502837,
              0.97189609, 0.28942311],
             [0.77261608, 0.6945811 , 0.30133295, ..., 0.04015439,
              0.0856194 , 0.33654333]],
    
            [[0.04720456, 0.11194704, 0.15198511, ..., 0.78657359,
              0.04491895, 0.97853686],
             [0.32533161, 0.42594342, 0.98040919, ..., 0.00241713,
              0.82341921, 0.3470506 ],
             [0.66521419, 0.91998776, 0.5643753 , ..., 0.65585365,
              0.62661315, 0.10590366],
             ...,
             [0.24925505, 0.05492279, 0.57887107, ..., 0.86056519,
              0.28726068, 0.82640873],
             [0.34793345, 0.12503262, 0.68898252, ..., 0.02282706,
              0.59485139, 0.73368563],
             [0.32162592, 0.21654125, 0.90026013, ..., 0.67946841,
              0.72849507, 0.2337043 ]]]], shape=(20, 20, 20, 20))

Quiz time again!

Can you guess the syntax?

quiz(2)

Getting the z slice is now easy and readable.

Adding coordinates

• Coordinates can be specified for each dimension.

• They describe the extent of each axis, as well as how far each data point is from its neighbours.

Here is an array that represents air pollution levels as a function of altitude and time.

data = sc.array(
    dims=["altitude", "year"],
    values=np.linspace(500, 10, 5).reshape((5, 1)) * rng.random(10),
)
sc.show(data)

data.plot()

In Scipp and Xarray, coordinates are added in a data structure called DataArray:

da = sc.DataArray(
    data=data,
    coords={
        "altitude": sc.linspace("altitude", 0, 8000, 5),
    },
)
sc.show(da)

da

scipp.DataArray (1.30 KB)

• Dimensions:
  • altitude: 5
  • year: 10
• Coordinates: (1)
  • altitude
    (altitude)
    float64
    𝟙
    0.0, 2000.000, 4000.000, 6000.000, 8000.000

    Values:
    array([   0., 2000., 4000., 6000., 8000.])

• Data: (1)
  • (altitude, year)
    float64
    𝟙
    488.350, 190.098, ..., 9.641, 2.636

    Values:
    array([[488.34988335, 190.09786751, 461.62311688, 130.84621193,
            159.54852921,  59.04561648, 120.88314663, 159.26696439,
            482.03962259, 131.82490214],
           [368.70416193, 143.52388997, 348.52545325,  98.78889001,
            120.45913955,  44.57944044,  91.2667757 , 120.24655812,
            363.93991505,  99.52780111],
           [249.05844051,  96.94991243, 235.42778961,  66.73156809,
             81.3697499 ,  30.11326441,  61.65040478,  81.22615184,
            245.84020752,  67.23070009],
           [129.41271909,  50.37593489, 122.33012597,  34.67424616,
             42.28036024,  15.64708837,  32.03403386,  42.20574556,
            127.74049999,  34.93359907],
           [  9.76699767,   3.80195735,   9.23246234,   2.61692424,
              3.19097058,   1.18091233,   2.41766293,   3.18533929,
              9.64079245,   2.63649804]])

da.plot()

Accessing and adding more coordinates

Coordinates are stored in a dict, and each dimension can have more than one coordinate.

Getting and setting coordinates is done using the same syntax as Python dicts:

print(da.coords.keys())
da.coords["altitude"]

<scipp.Dict.keys {altitude}>

scipp.Variable (296 Bytes)

• • (altitude: 5)
    float64
    𝟙
    0.0, 2000.000, 4000.000, 6000.000, 8000.000

    Values:
    array([   0., 2000., 4000., 6000., 8000.])

Exercise 1.1: Adding a new coordinate

The air pollution data was collected every year from 2014 to 2023; [2014, 2024).
Let’s add a coordinate, year to the year dimension.

Tip: You can create a Variable with consecutive numbers by using sc.arange(dim, start, stop).

Hint

da = sc.DataArray(
    data=data,
    coords={
        "altitude": sc.linspace("altitude", 0, 8000, 5),
        "year": sc.arange(..., 2014, ...)
    },
)

or

da.coords['year'] = sc.arange(..., 2014, ...)

Solution:

Show code cell content

Hide code cell content

da = sc.DataArray(
    data=data,
    coords={
        "altitude": sc.linspace("altitude", 0, 8000, 5),
        "year": sc.arange("year", 2014, 2024),
    },
)
sc.show(da)
da

scipp.DataArray (1.63 KB)

• Dimensions:
  • altitude: 5
  • year: 10
• Coordinates: (2)
  • altitude
    (altitude)
    float64
    𝟙
    0.0, 2000.000, 4000.000, 6000.000, 8000.000

    Values:
    array([   0., 2000., 4000., 6000., 8000.])

  • year
    (year)
    int64
    𝟙
    2014, 2015, ..., 2022, 2023

    Values:
    array([2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023])

• Data: (1)
  • (altitude, year)
    float64
    𝟙
    488.350, 190.098, ..., 9.641, 2.636

    Values:
    array([[488.34988335, 190.09786751, 461.62311688, 130.84621193,
            159.54852921,  59.04561648, 120.88314663, 159.26696439,
            482.03962259, 131.82490214],
           [368.70416193, 143.52388997, 348.52545325,  98.78889001,
            120.45913955,  44.57944044,  91.2667757 , 120.24655812,
            363.93991505,  99.52780111],
           [249.05844051,  96.94991243, 235.42778961,  66.73156809,
             81.3697499 ,  30.11326441,  61.65040478,  81.22615184,
            245.84020752,  67.23070009],
           [129.41271909,  50.37593489, 122.33012597,  34.67424616,
             42.28036024,  15.64708837,  32.03403386,  42.20574556,
            127.74049999,  34.93359907],
           [  9.76699767,   3.80195735,   9.23246234,   2.61692424,
              3.19097058,   1.18091233,   2.41766293,   3.18533929,
              9.64079245,   2.63649804]])

Exercise 1.2: Compute new coordinate

Add a new coordinate representing the Scipp-year.

Hint: Scipp was first released in 2020

Solution:

Show code cell content

Hide code cell content

da.coords["scipp-year"] = da.coords["year"] - 2020
sc.show(da)
da

scipp.DataArray (1.96 KB)

• Dimensions:
  • altitude: 5
  • year: 10
• Coordinates: (3)
  • altitude
    (altitude)
    float64
    𝟙
    0.0, 2000.000, 4000.000, 6000.000, 8000.000

    Values:
    array([   0., 2000., 4000., 6000., 8000.])

  • scipp-year
    (year)
    int64
    𝟙
    -6, -5, ..., 2, 3

    Values:
    array([-6, -5, -4, -3, -2, -1,  0,  1,  2,  3])

  • year
    (year)
    int64
    𝟙
    2014, 2015, ..., 2022, 2023

    Values:
    array([2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023])

• Data: (1)
  • (altitude, year)
    float64
    𝟙
    488.350, 190.098, ..., 9.641, 2.636

    Values:
    array([[488.34988335, 190.09786751, 461.62311688, 130.84621193,
            159.54852921,  59.04561648, 120.88314663, 159.26696439,
            482.03962259, 131.82490214],
           [368.70416193, 143.52388997, 348.52545325,  98.78889001,
            120.45913955,  44.57944044,  91.2667757 , 120.24655812,
            363.93991505,  99.52780111],
           [249.05844051,  96.94991243, 235.42778961,  66.73156809,
             81.3697499 ,  30.11326441,  61.65040478,  81.22615184,
            245.84020752,  67.23070009],
           [129.41271909,  50.37593489, 122.33012597,  34.67424616,
             42.28036024,  15.64708837,  32.03403386,  42.20574556,
            127.74049999,  34.93359907],
           [  9.76699767,   3.80195735,   9.23246234,   2.61692424,
              3.19097058,   1.18091233,   2.41766293,   3.18533929,
              9.64079245,   2.63649804]])

2. Going further

2.1 Physical units

Every data variable and coordinate in Scipp has physical units. (see also pint, astropy.units, pint-xarray)

Array Variable with unit:

temperature = sc.array(dims=["time"], values=[300.0, 301.0, 312.0, 340.0], unit="K")
temperature

scipp.Variable (288 Bytes)

• • (time: 4)
    float64
    K
    300.0, 301.0, 312.0, 340.0

    Values:
    array([300., 301., 312., 340.])

Scalar Variable (no dimensions) with unit:

sound_speed = sc.scalar(340.0, unit="m/s")
sound_speed

scipp.Variable (264 Bytes)

• • ()
    float64
    m/s
    340.0

    Values:
    array(340.)

Coordinates and data with units in a DataArray:

cph_air = sc.DataArray(
    data=sc.array(
        dims=["altitude", "year"],
        values=np.linspace(500, 10, 5).reshape((5, 1)) * rng.random(10),
        unit="m^-3",
    ),
    coords={
        "altitude": sc.linspace("altitude", 0, 8000, 5, unit="m"),
        "year": sc.arange("year", 2014, 2024, unit="year"),
    },
)
cph_air

scipp.DataArray (1.63 KB)

• Dimensions:
  • altitude: 5
  • year: 10
• Coordinates: (2)
  • altitude
    (altitude)
    float64
    m
    0.0, 2000.000, 4000.000, 6000.000, 8000.000

    Values:
    array([   0., 2000., 4000., 6000., 8000.])

  • year
    (year)
    int64
    Y
    2014, 2015, ..., 2022, 2023

    Values:
    array([2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023])

• Data: (1)
  • (altitude, year)
    float64
    1/m^3
    220.503, 304.935, ..., 1.721, 8.704

    Values:
    array([[220.50306103, 304.93540471, 431.81064828, 431.87883539,
            337.44065667, 329.93717398, 367.87884916, 111.37682907,
             86.03309233, 435.20748624],
           [166.47981108, 230.22623056, 326.01703946, 326.06852072,
            254.76769579, 249.10256635, 277.74853111,  84.08950594,
             64.95498471, 328.58165211],
           [112.45656112, 155.5170564 , 220.22343063, 220.25820605,
            172.0947349 , 168.26795873, 187.61821307,  56.80218282,
             43.87687709, 221.95581798],
           [ 58.43331117,  80.80788225, 114.4298218 , 114.44789138,
             89.42177402,  87.4333511 ,  97.48789503,  29.5148597 ,
             22.79876947, 115.32998385],
           [  4.41006122,   6.09870809,   8.63621297,   8.63757671,
              6.74881313,   6.59874348,   7.35757698,   2.22753658,
              1.72066185,   8.70414972]])

Units are automatically handled in arithmetic operations.

Say we know the mean ultra-fine particle mass

ultra_fine_particle_mass = sc.scalar(1.0e-6, unit="kg")

cph_air *= ultra_fine_particle_mass
cph_air

scipp.DataArray (1.63 KB)

• Dimensions:
  • altitude: 5
  • year: 10
• Coordinates: (2)
  • altitude
    (altitude)
    float64
    m
    0.0, 2000.000, 4000.000, 6000.000, 8000.000

    Values:
    array([   0., 2000., 4000., 6000., 8000.])

  • year
    (year)
    int64
    Y
    2014, 2015, ..., 2022, 2023

    Values:
    array([2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023])

• Data: (1)
  • (altitude, year)
    float64
    kg/m^3
    0.000, 0.000, ..., 1.721e-06, 8.704e-06

    Values:
    array([[2.20503061e-04, 3.04935405e-04, 4.31810648e-04, 4.31878835e-04,
            3.37440657e-04, 3.29937174e-04, 3.67878849e-04, 1.11376829e-04,
            8.60330923e-05, 4.35207486e-04],
           [1.66479811e-04, 2.30226231e-04, 3.26017039e-04, 3.26068521e-04,
            2.54767696e-04, 2.49102566e-04, 2.77748531e-04, 8.40895059e-05,
            6.49549847e-05, 3.28581652e-04],
           [1.12456561e-04, 1.55517056e-04, 2.20223431e-04, 2.20258206e-04,
            1.72094735e-04, 1.68267959e-04, 1.87618213e-04, 5.68021828e-05,
            4.38768771e-05, 2.21955818e-04],
           [5.84333112e-05, 8.08078822e-05, 1.14429822e-04, 1.14447891e-04,
            8.94217740e-05, 8.74333511e-05, 9.74878950e-05, 2.95148597e-05,
            2.27987695e-05, 1.15329984e-04],
           [4.41006122e-06, 6.09870809e-06, 8.63621297e-06, 8.63757671e-06,
            6.74881313e-06, 6.59874348e-06, 7.35757698e-06, 2.22753658e-06,
            1.72066185e-06, 8.70414972e-06]])

Units also provide protection

Say we now also have air pollution data for another city, e.g., NYC.

We would like to compute the difference between CPH and NYC air pollution (as a function of altitude and year), but we forgot to multiply the NYC data by particle mass:

nyc_air = sc.DataArray(
    data=sc.array(
        dims=["altitude", "year"],
        values=np.linspace(800, 20, 5).reshape((5, 1)) * rng.random(10),
        unit="m-3",
    ),
    coords={
        "altitude": sc.linspace("altitude", 0, 8000, 5, unit="m"),
        "year": sc.arange("year", 2014, 2024, unit="year"),
    },
)

cph_air - nyc_air

---------------------------------------------------------------------------
UnitError                                 Traceback (most recent call last)
Cell In[24], line 13
      9         "year": sc.arange("year", 2014, 2024, unit="year"),
     10     },
     11 )
     12 
---> 13 cph_air - nyc_air

UnitError: Cannot subtract kg/m^3 and 1/m^3.

nyc_air *= ultra_fine_particle_mass

air_difference = cph_air - nyc_air

air_difference.plot()

• Units are very useful in early prevention of difficult-to-spot bugs in a workflow.

• They save hours of debugging time, free-up mental capacity and let the user focus on the important thing: doing science.

Units for label-based indexing

We also use units to distinguish between positional indexing and label-based indexing:

cph_air["altitude", 2000.0 * sc.Unit("m")].plot()

Positional indices are based on the dimension, and value indices are based on the coordinates.

Exercise 2: Coordinate and Units

We have a data array that contains air pollution as a function of year and altitude above the city of Copenhagen. However, we want to have a pressure coordinate for the altitude dimension instead of altitude.

Assuming a constant air temperature \(T\) of 300 K, the pressure as a function of height \(h\) is given by

\[ P = P_{0} \exp{\left[ \frac{-g_{0}Mh}{RT} \right]} \]

Here is the incomplete function altitude_to_pressure that converts altitude[m] into pressure[hPa].

Complete the function and use it to add the pressure coordinate to cph_air.

def altitude_to_pressure(altitude):
    M = sc.scalar(0.0289644, unit="kg/mol")
    g0 = sc.scalar(9.80665, unit="m/s2")
    R = sc.scalar(8.3144598, unit="J/mol/K")
    T = sc.scalar(300.0)
    p0 = sc.scalar(1013.25, unit="hPa")
    return p0 * sc.exp(-g0 * M * altitude / (R * T))

Solution:

Show code cell content

Hide code cell content

def altitude_to_pressure(altitude):
    M = sc.scalar(0.0289644, unit="kg/mol")
    g0 = sc.scalar(9.80665, unit="m/s2")
    R = sc.scalar(8.3144598, unit="J/mol/K")
    T = sc.scalar(300.0, unit="K")
    p0 = sc.scalar(1013.25, unit="hPa")
    return p0 * sc.exp(-g0 * M * altitude / (R * T))


cph_air.coords["pressure"] = altitude_to_pressure(cph_air.coords["altitude"])
cph_air

scipp.DataArray (1.92 KB)

• Dimensions:
  • altitude: 5
  • year: 10
• Coordinates: (3)
  • altitude
    (altitude)
    float64
    m
    0.0, 2000.000, 4000.000, 6000.000, 8000.000

    Values:
    array([   0., 2000., 4000., 6000., 8000.])

  • pressure
    (altitude)
    float64
    100Pa
    1013.250, 806.874, 642.532, 511.663, 407.449

    Values:
    array([1013.25      ,  806.87395253,  642.53202593,  511.66282298,
            407.44870895])

  • year
    (year)
    int64
    Y
    2014, 2015, ..., 2022, 2023

    Values:
    array([2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023])

• Data: (1)
  • (altitude, year)
    float64
    kg/m^3
    0.000, 0.000, ..., 1.721e-06, 8.704e-06

    Values:
    array([[2.20503061e-04, 3.04935405e-04, 4.31810648e-04, 4.31878835e-04,
            3.37440657e-04, 3.29937174e-04, 3.67878849e-04, 1.11376829e-04,
            8.60330923e-05, 4.35207486e-04],
           [1.66479811e-04, 2.30226231e-04, 3.26017039e-04, 3.26068521e-04,
            2.54767696e-04, 2.49102566e-04, 2.77748531e-04, 8.40895059e-05,
            6.49549847e-05, 3.28581652e-04],
           [1.12456561e-04, 1.55517056e-04, 2.20223431e-04, 2.20258206e-04,
            1.72094735e-04, 1.68267959e-04, 1.87618213e-04, 5.68021828e-05,
            4.38768771e-05, 2.21955818e-04],
           [5.84333112e-05, 8.08078822e-05, 1.14429822e-04, 1.14447891e-04,
            8.94217740e-05, 8.74333511e-05, 9.74878950e-05, 2.95148597e-05,
            2.27987695e-05, 1.15329984e-04],
           [4.41006122e-06, 6.09870809e-06, 8.63621297e-06, 8.63757671e-06,
            6.74881313e-06, 6.59874348e-06, 7.35757698e-06, 2.22753658e-06,
            1.72066185e-06, 8.70414972e-06]])

2.2 Histogramming and bin-edge coordinates

• It is sometimes necessary to have coordinates that represent a range for each data value.

• E.g., “the temperature was 310 K in the time span between 10 and 20 seconds”.

• This also arises every time we histogram data.

• Scipp supports this by having bin-edge coordinates: a coordinate which has a length of 1 more than the dimension length.

The next data set is meant to represent photon events in a camera. We have a long list of x and y positions for the photons.

x = sc.array(dims=["row"], values=rng.normal(size=10000), unit="cm")
y = sc.array(dims=["row"], values=rng.normal(size=10000), unit="cm")
recording = sc.DataArray(
    data=sc.ones(sizes=x.sizes, unit="counts"), coords={"x": x, "y": y}
)
recording

scipp.DataArray (235.49 KB)

• Dimensions:
  • row: 10000
• Coordinates: (2)
  • x
    (row)
    float64
    cm
    -0.679, -0.621, ..., -0.788, 1.123

    Values:
    array([-0.67924997, -0.62053203,  1.33121422, ...,  0.66931119,
           -0.78846271,  1.12268421], shape=(10000,))

  • y
    (row)
    float64
    cm
    -0.306, -0.704, ..., -0.848, 0.739

    Values:
    array([-0.30610691, -0.70374693, -1.00885753, ..., -0.42935023,
           -0.84826518,  0.73921398], shape=(10000,))

• Data: (1)
  • (row)
    float64
    counts
    1.0, 1.0, ..., 1.0, 1.0

    Values:
    array([1., 1., 1., ..., 1., 1., 1.], shape=(10000,))

scatter(x.values, y.values)

It is very common to histogram such data.

In Scipp, histogramming has a very concise and easy-to-use syntax. To make 8 bins in both the x and y dimensions:

image = recording.hist(y=8, x=8)
image.plot(aspect="equal")

The x and y coordinates are now bin-edge coordinates.

sc.show(image)
image

scipp.DataArray (1.76 KB)

• Dimensions:
  • y: 8
  • x: 8
• Coordinates: (2)
  • x
    (x [bin-edge])
    float64
    cm
    -3.819, -2.831, ..., 3.095, 4.083

    Values:
    array([-3.81886182, -2.83110929, -1.84335675, -0.85560421,  0.13214833,
            1.11990087,  2.10765341,  3.09540594,  4.08315848])

  • y
    (y [bin-edge])
    float64
    cm
    -3.646, -2.663, ..., 3.236, 4.219

    Values:
    array([-3.64602432, -2.66291545, -1.67980659, -0.69669773,  0.28641113,
            1.26952   ,  2.25262886,  3.23573772,  4.21884658])

• Data: (1)
  • (y, x)
    float64
    counts
    0.0, 2.0, ..., 1.0, 0.0

    Values:
    array([[0.000e+00, 2.000e+00, 8.000e+00, 1.200e+01, 1.000e+01, 5.000e+00,
            0.000e+00, 0.000e+00],
           [1.000e+00, 1.100e+01, 8.100e+01, 1.460e+02, 1.190e+02, 4.600e+01,
            6.000e+00, 1.000e+00],
           [4.000e+00, 6.200e+01, 3.290e+02, 7.620e+02, 6.310e+02, 2.370e+02,
            3.000e+01, 3.000e+00],
           [7.000e+00, 1.210e+02, 5.790e+02, 1.298e+03, 1.179e+03, 4.170e+02,
            7.000e+01, 2.000e+00],
           [5.000e+00, 8.800e+01, 4.440e+02, 1.030e+03, 8.770e+02, 3.240e+02,
            4.400e+01, 1.000e+00],
           [3.000e+00, 3.600e+01, 1.490e+02, 3.280e+02, 2.710e+02, 9.300e+01,
            1.400e+01, 2.000e+00],
           [0.000e+00, 5.000e+00, 1.900e+01, 3.500e+01, 3.600e+01, 1.100e+01,
            1.000e+00, 0.000e+00],
           [0.000e+00, 0.000e+00, 1.000e+00, 2.000e+00, 1.000e+00, 0.000e+00,
            1.000e+00, 0.000e+00]])

• Numpy and Matplotlib return the bin edges and the data counts separately.

• We have everything stored inside a single data structure.

You can, of course, adjust the number of bins:

recording.hist(y=100, x=100).plot(aspect="equal")

Exercise 3: Histogramming

We found a 2D detector that reads your mood!

We recorded a signal with it, and now we can visualize the signal by histogramming.

from scipp_utils import load_signal_to_histogram

signal_rng = np.random.default_rng(1)
signal = load_signal_to_histogram(signal_rng)
signal

scipp.DataArray (4.66 MB)

• Dimensions:
  • row: 203700
• Coordinates: (2)
  • x
    (row)
    float64
    cm
    7.278, -16.296, ..., 0.018, 0.016

    Values:
    array([ 7.27805283e+00, -1.62959448e+01, -7.09559402e+00, ...,
            1.89787051e-02,  1.81696472e-02,  1.60307256e-02], shape=(203700,))

  • y
    (row)
    float64
    cm
    0.432, -0.784, ..., 0.061, 0.060

    Values:
    array([ 0.43193503, -0.78422306, -9.62131874, ...,  0.06044533,
            0.06091188,  0.05965346], shape=(203700,))

• Data: (1)
  • (row)
    float64
    counts
    1.0, 1.0, ..., 1.0, 1.0

    Values:
    array([1., 1., 1., ..., 1., 1., 1.], shape=(203700,))

Exercise 3-1: Number of bins for histogramming.

First, we need to find the right number of bins to histogram the signal.

We tried 200 bins and 4 bins for each axis, but none of them seems meaningful!

signal.hist(x=200, y=200).plot() + signal.hist(x=4, y=4).plot()

Solution:

Show code cell content

Hide code cell content

# 30~50 bins are enough to see the meaningful shape!
signal.hist(x=50, y=50).plot() + signal.hist(x=30, y=30).plot()

Exercise 3-2: Custom histogram edges.

However, there is a suspicious hot spot in the very middle of the image.

We want to investigate those signals within the specific range of x and y.

Let’s histogram the hot spot and see what is in there.

You can histogram the data with custom histogram edges like below.

Hint:

hist_edges_x = sc.linspace(dim="x", start=-10, stop=10, unit="cm", num=200)
hist_edges_y = sc.linspace(dim="y", start=-10, stop=10, unit="cm", num=200)
signal.hist(x=hist_edges_x, y=hist_edges_y).plot()

Solution:

Show code cell content

Hide code cell content

# There was a smiley in the middle of the heart!

hist_edges_x = sc.linspace(dim="x", start=-0.15, stop=0.15, unit="cm", num=50)
hist_edges_y = sc.linspace(dim="y", start=-0.15, stop=0.15, unit="cm", num=50)
signal.hist(x=hist_edges_x, y=hist_edges_y).plot()

3. Binned data

Scipp distinguishes histogrammed data from binned data:

• Histogrammed data refers to regular dense arrays of, e.g., floating-point values with an associated bin-edge coordinate.

• Binned data refers to the precursor of histogrammed data, i.e., each bin contains a “list” of contributing events or values. Binned data can be converted into a histogram by computing the sum over all events or values in a bin.

This is conceptually similar to a multi-dimensional .

It is best illustrated with an example of data analysis. For this, we will use a dataset containing simulated neutron events.

A special sample for imaging

file = fetch_data("4-reduction/special_imaging_sample_small")

da = sc.io.load_hdf5(file)
da

scipp.DataArray (1.48 GB)

• Dimensions:
  • event: 13240308
• Coordinates: (14)
  • birth_time
    (event)
    float64
    µs
    2095.575, 1629.573, ..., 1330.540, 994.222

    Values:
    array([2095.57544851, 1629.57329071, 1415.45174112, ..., 2491.22896273,
           1330.53994267,  994.22156859], shape=(13240308,))

  • birth_x
    (event)
    float64
    m
    -0.052, -0.103, ..., 0.073, -0.089

    Values:
    array([-0.05158719, -0.10339723,  0.07845079, ..., -0.0560222 ,
            0.07335051, -0.0888617 ], shape=(13240308,))

  • birth_y
    (event)
    float64
    m
    0.029, -0.009, ..., 0.013, -0.013

    Values:
    array([ 0.02896936, -0.00869429,  0.00833677, ..., -0.00641298,
            0.01259691, -0.0133081 ], shape=(13240308,))

  • birth_z
    (event)
    float64
    m
    0.021, 0.013, ..., 0.021, 0.015

    Values:
    array([0.02054448, 0.01275013, 0.01774066, ..., 0.02276014, 0.02068504,
           0.01510914], shape=(13240308,))

  • energy_transfer
    (event)
    float64
    meV
    0.0, 0.0, ..., 0.0, 0.0

    Values:
    array([0., 0., 0., ..., 0., 0., 0.], shape=(13240308,))

  • final_energy
    (event)
    float64
    meV
    2.648, 0.632, ..., 2.743, 0.375

    Values:
    array([ 2.64844075,  0.63159722, 10.85213099, ...,  0.66819168,
            2.74277137,  0.3752962 ], shape=(13240308,))

  • initial_energy
    (event)
    float64
    meV
    2.648, 0.632, ..., 2.743, 0.375

    Values:
    array([ 2.64844075,  0.63159722, 10.85213099, ...,  0.66819168,
            2.74277137,  0.3752962 ], shape=(13240308,))

  • initial_wavelength
    (event)
    float64
    Å
    5.558, 11.381, ..., 5.461, 14.764

    Values:
    array([ 5.55766926, 11.3806672 ,  2.74555605, ..., 11.06464022,
            5.46126216, 14.76388809], shape=(13240308,))

  • reflections
    (event)
    int64
    𝟙
    3, 13, ..., 0, 18

    Values:
    array([ 3, 13,  1, ...,  8,  0, 18], shape=(13240308,))

  • toa
    (event)
    float64
    µs
    8.709e+04, 1.757e+05, ..., 8.485e+04, 2.268e+05

    Values:
    array([ 87089.53380614, 175675.1865625 ,  43403.49855428, ...,
           171703.81184739,  84850.13550796, 226779.74993476],
          shape=(13240308,))

  • wavelength
    (event)
    float64
    Å
    5.558, 11.381, ..., 5.461, 14.764

    Values:
    array([ 5.55766926, 11.3806672 ,  2.74555605, ..., 11.06464022,
            5.46126216, 14.76388809], shape=(13240308,))

  • x
    (event)
    float64
    m
    -0.754, -0.812, ..., 0.775, 0.606

    Values:
    array([-0.7542735 , -0.81190604, -0.80441472, ...,  0.71049372,
            0.77481493,  0.6058382 ], shape=(13240308,))

  • y
    (event)
    float64
    m
    -0.701, -0.671, ..., 0.700, 0.719

    Values:
    array([-0.70125091, -0.67064033, -0.63500747, ...,  0.84913847,
            0.70014981,  0.71859393], shape=(13240308,))

  • z
    (event)
    float64
    m
    60.5, 60.5, ..., 60.5, 60.5

    Values:
    array([60.5, 60.5, 60.5, ..., 60.5, 60.5, 60.5], shape=(13240308,))

• Data: (1)
  • (event)
    float64
    𝟙
    1.0, 1.0, ..., 1.0, 1.0

    Values:
    array([1., 1., 1., ..., 1., 1., 1.], shape=(13240308,))

n = 50
x = da.coords["x"].values[::n]
y = da.coords["y"].values[::n]
scatter(x, y)

Binning the data records

• Working with binned data is most efficient when keeping the number of bins relatively low.

• Binning is essentially like overlaying a grid of bin edges onto our data

ax = scatter(x, y, get_ax=True)
for xline in np.linspace(*ax.get_xlim(), 7):
    ax.axvline(xline, color="k")
for yline in np.linspace(*ax.get_ylim(), 7):
    ax.axhline(yline, color="k")

# Bin into 6 x & y bins
binned = da.bin(y=6, x=6)
binned

scipp.DataArray (1.48 GB)

• Dimensions:
  • y: 6
  • x: 6
• Coordinates: (2)
  • x
    (x [bin-edge])
    float64
    m
    -0.888, -0.592, ..., 0.592, 0.888

    Values:
    array([-8.87999982e-01, -5.91999991e-01, -2.96000000e-01, -9.49421353e-09,
            2.95999981e-01,  5.91999972e-01,  8.87999963e-01])

  • y
    (y [bin-edge])
    float64
    m
    -0.888, -0.592, ..., 0.592, 0.888

    Values:
    array([-8.87999565e-01, -5.91999802e-01, -2.96000038e-01, -2.75037134e-07,
            2.95999488e-01,  5.91999252e-01,  8.87999015e-01])

• Data: (1)
  • (y, x)
    float64
    𝟙
    binned data [len=415497, len=403963, ..., len=395036, len=234432]

    dim='event',
    content=DataArray(
              dims=(event: 13240308),
              data=float64[dimensionless],
              coords={'x':float64[m], 'y':float64[m], 'z':float64[m], 'wavelength':float64[Å],
                      'initial_energy':float64[meV], 'initial_wavelength':float64[Å],
                      'final_energy':float64[meV], 'energy_transfer':float64[meV],
                      'reflections':int64[dimensionless], 'birth_x':float64[m],
                      'birth_y':float64[m], 'birth_z':float64[m], 'birth_time':float64[µs],
                      'toa':float64[µs]})

# Histogramming is summing all the counts in each bin
binned_sum = binned.bins.sum()

binned_sum.plot(aspect="equal")

Selecting/slicing bins

• Binning groups the data into bins, but keeps the underlying table of records.

• No information is lost, it is simply re-ordered.

• The bins can then be used for slicing the data, providing extremely efficient data selection and filtering.

mat_a = binned["x", 1]["y", 5]
mat_a

scipp.DataArray (42.45 MB out of 1.48 GB)

• Dimensions:
• Coordinates: (2)
  • x
    (x [bin-edge])
    float64
    m
    -0.592, -0.296

    Values:
    array([-0.59199999, -0.296     ])

  • y
    (y [bin-edge])
    float64
    m
    0.592, 0.888

    Values:
    array([0.59199925, 0.88799901])

• Data: (1)
  • ()
    float64
    𝟙
    binned data [len=370874]

    dim='event',
    content=DataArray(
              dims=(event: 13240308),
              data=float64[dimensionless],
              coords={'x':float64[m], 'y':float64[m], 'z':float64[m], 'wavelength':float64[Å],
                      'initial_energy':float64[meV], 'initial_wavelength':float64[Å],
                      'final_energy':float64[meV], 'energy_transfer':float64[meV],
                      'reflections':int64[dimensionless], 'birth_x':float64[m],
                      'birth_y':float64[m], 'birth_z':float64[m], 'birth_time':float64[µs],
                      'toa':float64[µs]})

# We can now histogram this with a much finer resolution

mat_a.hist(y=300, x=300).plot(aspect="equal")

# We select another bin, which contains a second material

mat_b = binned["x", 0]["y", 1]
mat_b.hist(y=300, x=300).plot(aspect="equal")

Binning into a new dimension: wavelength spectra

• Data that has already been binned can also be binned further into new dimensions.

Fe = mat_a
Al = mat_b

Fe

scipp.DataArray (42.45 MB out of 1.48 GB)

• Dimensions:
• Coordinates: (2)
  • x
    (x [bin-edge])
    float64
    m
    -0.592, -0.296

    Values:
    array([-0.59199999, -0.296     ])

  • y
    (y [bin-edge])
    float64
    m
    0.592, 0.888

    Values:
    array([0.59199925, 0.88799901])

• Data: (1)
  • ()
    float64
    𝟙
    binned data [len=370874]

    dim='event',
    content=DataArray(
              dims=(event: 13240308),
              data=float64[dimensionless],
              coords={'x':float64[m], 'y':float64[m], 'z':float64[m], 'wavelength':float64[Å],
                      'initial_energy':float64[meV], 'initial_wavelength':float64[Å],
                      'final_energy':float64[meV], 'energy_transfer':float64[meV],
                      'reflections':int64[dimensionless], 'birth_x':float64[m],
                      'birth_y':float64[m], 'birth_z':float64[m], 'birth_time':float64[µs],
                      'toa':float64[µs]})

• We look at the wavelength spectrum inside the Fe and Al bins we have selected above.

# Use 300 wavelength bins
Fe_spectrum = Fe.bin(wavelength=300)
Fe_spectrum

scipp.DataArray (42.46 MB out of 42.46 MB)

• Dimensions:
  • wavelength: 300
• Coordinates: (3)
  • wavelength
    (wavelength [bin-edge])
    float64
    Å
    0.050, 0.100, ..., 14.950, 15.000

    Values:
    array([ 0.05006454,  0.09989759,  0.14973064,  0.19956369,  0.24939675,
            0.2992298 ,  0.34906285,  0.3988959 ,  0.44872895,  0.498562  ,
            0.54839505,  0.59822811,  0.64806116,  0.69789421,  0.74772726,
            0.79756031,  0.84739336,  0.89722642,  0.94705947,  0.99689252,
            1.04672557,  1.09655862,  1.14639167,  1.19622472,  1.24605778,
            1.29589083,  1.34572388,  1.39555693,  1.44538998,  1.49522303,
            1.54505608,  1.59488914,  1.64472219,  1.69455524,  1.74438829,
            1.79422134,  1.84405439,  1.89388744,  1.9437205 ,  1.99355355,
            2.0433866 ,  2.09321965,  2.1430527 ,  2.19288575,  2.24271881,
            2.29255186,  2.34238491,  2.39221796,  2.44205101,  2.49188406,
            2.54171711,  2.59155017,  2.64138322,  2.69121627,  2.74104932,
            2.79088237,  2.84071542,  2.89054847,  2.94038153,  2.99021458,
            3.04004763,  3.08988068,  3.13971373,  3.18954678,  3.23937984,
            3.28921289,  3.33904594,  3.38887899,  3.43871204,  3.48854509,
            3.53837814,  3.5882112 ,  3.63804425,  3.6878773 ,  3.73771035,
            3.7875434 ,  3.83737645,  3.8872095 ,  3.93704256,  3.98687561,
            4.03670866,  4.08654171,  4.13637476,  4.18620781,  4.23604086,
            4.28587392,  4.33570697,  4.38554002,  4.43537307,  4.48520612,
            4.53503917,  4.58487223,  4.63470528,  4.68453833,  4.73437138,
            4.78420443,  4.83403748,  4.88387053,  4.93370359,  4.98353664,
            5.03336969,  5.08320274,  5.13303579,  5.18286884,  5.23270189,
            5.28253495,  5.332368  ,  5.38220105,  5.4320341 ,  5.48186715,
            5.5317002 ,  5.58153326,  5.63136631,  5.68119936,  5.73103241,
            5.78086546,  5.83069851,  5.88053156,  5.93036462,  5.98019767,
            6.03003072,  6.07986377,  6.12969682,  6.17952987,  6.22936292,
            6.27919598,  6.32902903,  6.37886208,  6.42869513,  6.47852818,
            6.52836123,  6.57819429,  6.62802734,  6.67786039,  6.72769344,
            6.77752649,  6.82735954,  6.87719259,  6.92702565,  6.9768587 ,
            7.02669175,  7.0765248 ,  7.12635785,  7.1761909 ,  7.22602395,
            7.27585701,  7.32569006,  7.37552311,  7.42535616,  7.47518921,
            7.52502226,  7.57485531,  7.62468837,  7.67452142,  7.72435447,
            7.77418752,  7.82402057,  7.87385362,  7.92368668,  7.97351973,
            8.02335278,  8.07318583,  8.12301888,  8.17285193,  8.22268498,
            8.27251804,  8.32235109,  8.37218414,  8.42201719,  8.47185024,
            8.52168329,  8.57151634,  8.6213494 ,  8.67118245,  8.7210155 ,
            8.77084855,  8.8206816 ,  8.87051465,  8.92034771,  8.97018076,
            9.02001381,  9.06984686,  9.11967991,  9.16951296,  9.21934601,
            9.26917907,  9.31901212,  9.36884517,  9.41867822,  9.46851127,
            9.51834432,  9.56817737,  9.61801043,  9.66784348,  9.71767653,
            9.76750958,  9.81734263,  9.86717568,  9.91700873,  9.96684179,
           10.01667484, 10.06650789, 10.11634094, 10.16617399, 10.21600704,
           10.2658401 , 10.31567315, 10.3655062 , 10.41533925, 10.4651723 ,
           10.51500535, 10.5648384 , 10.61467146, 10.66450451, 10.71433756,
           10.76417061, 10.81400366, 10.86383671, 10.91366976, 10.96350282,
           11.01333587, 11.06316892, 11.11300197, 11.16283502, 11.21266807,
           11.26250113, 11.31233418, 11.36216723, 11.41200028, 11.46183333,
           11.51166638, 11.56149943, 11.61133249, 11.66116554, 11.71099859,
           11.76083164, 11.81066469, 11.86049774, 11.91033079, 11.96016385,
           12.0099969 , 12.05982995, 12.109663  , 12.15949605, 12.2093291 ,
           12.25916216, 12.30899521, 12.35882826, 12.40866131, 12.45849436,
           12.50832741, 12.55816046, 12.60799352, 12.65782657, 12.70765962,
           12.75749267, 12.80732572, 12.85715877, 12.90699182, 12.95682488,
           13.00665793, 13.05649098, 13.10632403, 13.15615708, 13.20599013,
           13.25582318, 13.30565624, 13.35548929, 13.40532234, 13.45515539,
           13.50498844, 13.55482149, 13.60465455, 13.6544876 , 13.70432065,
           13.7541537 , 13.80398675, 13.8538198 , 13.90365285, 13.95348591,
           14.00331896, 14.05315201, 14.10298506, 14.15281811, 14.20265116,
           14.25248421, 14.30231727, 14.35215032, 14.40198337, 14.45181642,
           14.50164947, 14.55148252, 14.60131558, 14.65114863, 14.70098168,
           14.75081473, 14.80064778, 14.85048083, 14.90031388, 14.95014694,
           14.99997999])

  • x
    (x [bin-edge])
    float64
    m
    -0.592, -0.296

    Values:
    array([-0.59199999, -0.296     ])

  • y
    (y [bin-edge])
    float64
    m
    0.592, 0.888

    Values:
    array([0.59199925, 0.88799901])

• Data: (1)
  • (wavelength)
    float64
    𝟙
    binned data [len=832, len=810, ..., len=1347, len=1358]

    dim='event',
    content=DataArray(
              dims=(event: 370874),
              data=float64[dimensionless],
              coords={'x':float64[m], 'y':float64[m], 'z':float64[m], 'wavelength':float64[Å],
                      'initial_energy':float64[meV], 'initial_wavelength':float64[Å],
                      'final_energy':float64[meV], 'energy_transfer':float64[meV],
                      'reflections':int64[dimensionless], 'birth_x':float64[m],
                      'birth_y':float64[m], 'birth_z':float64[m], 'birth_time':float64[µs],
                      'toa':float64[µs]})

Fe_spectrum.hist().plot()

Al_spectrum = Al.bin(wavelength=300)

# Plot both spectra on the same axes
pp.plot({"Fe": Fe_spectrum.hist(), "Al": Al_spectrum.hist()})

We can also compare this to the same plot for the entire sample slab:

da.hist(wavelength=300).plot()

Other operations on bins: are there any materials that exhibit inelastic behaviour?

• In addition to summing/histogramming, bins can be used for other reduction operations: min(), max(), and mean().

binned

scipp.DataArray (1.48 GB)

• Dimensions:
  • y: 6
  • x: 6
• Coordinates: (2)
  • x
    (x [bin-edge])
    float64
    m
    -0.888, -0.592, ..., 0.592, 0.888

    Values:
    array([-8.87999982e-01, -5.91999991e-01, -2.96000000e-01, -9.49421353e-09,
            2.95999981e-01,  5.91999972e-01,  8.87999963e-01])

  • y
    (y [bin-edge])
    float64
    m
    -0.888, -0.592, ..., 0.592, 0.888

    Values:
    array([-8.87999565e-01, -5.91999802e-01, -2.96000038e-01, -2.75037134e-07,
            2.95999488e-01,  5.91999252e-01,  8.87999015e-01])

• Data: (1)
  • (y, x)
    float64
    𝟙
    binned data [len=415497, len=403963, ..., len=395036, len=234432]

    dim='event',
    content=DataArray(
              dims=(event: 13240308),
              data=float64[dimensionless],
              coords={'x':float64[m], 'y':float64[m], 'z':float64[m], 'wavelength':float64[Å],
                      'initial_energy':float64[meV], 'initial_wavelength':float64[Å],
                      'final_energy':float64[meV], 'energy_transfer':float64[meV],
                      'reflections':int64[dimensionless], 'birth_x':float64[m],
                      'birth_y':float64[m], 'birth_z':float64[m], 'birth_time':float64[µs],
                      'toa':float64[µs]})

• To get the maximum energy transfer for each material, we simply compute the max() of energy transfer inside each bin:

binned.bins.coords["energy_transfer"].bins.max()

scipp.Variable (544 Bytes)

• • (y: 6, x: 6)
    float64
    meV
    0.0, 0.0, ..., 0.0, 0.0

    Values:
    array([[0.        , 0.        , 0.        , 0.        , 0.        ,
            0.        ],
           [0.        , 0.        , 0.        , 0.        , 0.        ,
            0.        ],
           [0.1219953 , 0.        , 0.        , 0.        , 0.        ,
            0.        ],
           [0.        , 0.        , 0.        , 0.        , 0.06286986,
            0.        ],
           [0.017002  , 0.        , 0.06565547, 0.        , 0.        ,
            0.        ],
           [0.        , 0.        , 0.        , 0.        , 0.        ,
            0.        ]])

• Using the values, we can identify the bins that show evidence of inelastic processes.

binned["y", 2]["x", 0].hist(energy_transfer=300).plot()

binned["y", 4]["x", 0].hist(energy_transfer=300).plot()

Reflections

Another interesting quantity to look at is the number of reflections along the guide as a function of wavelength. We first inspect this for all the events in the data array:

fig = da.hist(reflections=200, wavelength=300).plot(logc=True, cmap="magma")
fig

This is telling us that lower energy neutrons are able to do more reflections: this is because the critical angle for total reflection is proportional to wavelength, so the guide accepts more incident angle trajectories for longer wavelengths.

We then select one bin at the edge of the sample, and one bin close to the center, and overlay them as scatter plots:

edge = binned["y", 0]["x", 0].value
center = binned["y", 3]["x", 3].value

pp.scatter(
    {"edge": edge[::100], "center": center[::100]},
    x="wavelength",
    y="reflections",
    marker=".",
    color={"edge": "lime", "center": "C0"},
    ax=fig.ax,
)

This reveals that neutrons hitting the edge of the sample have undergone slightly more reflections than the ones arriving in the middle.

4. More data visualizations with Plopp

https://scipp.github.io/plopp

The slicer plot

hist_e = da.hist(energy_transfer=256, y=256, x=256)

pp.slicer(hist_e, aspect="equal", autoscale=False)

The inspector plot

hist_w = da.hist(wavelength=256, y=256, x=256)

pp.inspector(
    hist_w, dim="wavelength", aspect="equal", autoscale=False, mode="rectangle"
)

The scatter3d plot

Looking at the birth positions of the neutrons, we reveal the surface of the ESS moderator, and wavelength anisotropies.

subset = da[::100].copy()
subset.data = subset.coords["wavelength"]

pp.scatter3d(subset, x="birth_x", y="birth_y", z="birth_z", size=0.001, cbar=True)