12. Creating a DataStore from numpy arrays

The goal of this notebook is to demonstrate how to create a DataStore from scratch. This can be useful if your device is not supported or if you would like to integrate the dtscalibration library in your current routine.

import numpy as np
import os
import matplotlib.pyplot as plt
import xarray as xr

from dtscalibration import DataStore, read_silixa_files

For a DataStore object, a few things are needed:

• timestamps

• Stokes signal

• anti-Stokes signal

• x (length along fiber)

Let’s grab the data from an existing silixa dataset:

filepath = os.path.join("..", "..", "tests", "data", "single_ended")

ds_silixa = read_silixa_files(directory=filepath, silent=True)

/home/docs/checkouts/readthedocs.org/user_builds/python-dts-calibration/envs/latest/lib/python3.9/site-packages/dtscalibration/datastore.py:81: UserWarning: Converting non-nanosecond precision timedelta values to nanosecond precision. This behavior can eventually be relaxed in xarray, as it is an artifact from pandas which is now beginning to support non-nanosecond precision values. This warning is caused by passing non-nanosecond np.datetime64 or np.timedelta64 values to the DataArray or Variable constructor; it can be silenced by converting the values to nanosecond precision ahead of time.
  super().__init__(*args, **kwargs)
/home/docs/checkouts/readthedocs.org/user_builds/python-dts-calibration/envs/latest/lib/python3.9/site-packages/dtscalibration/datastore.py:81: UserWarning: Converting non-nanosecond precision timedelta values to nanosecond precision. This behavior can eventually be relaxed in xarray, as it is an artifact from pandas which is now beginning to support non-nanosecond precision values. This warning is caused by passing non-nanosecond np.datetime64 or np.timedelta64 values to the DataArray or Variable constructor; it can be silenced by converting the values to nanosecond precision ahead of time.
  super().__init__(*args, **kwargs)

We will get all the numpy arrays from this DataStore to create a new one from ‘scratch’.

Let’s start with the most basic data:

x = ds_silixa.x.values
time = ds_silixa.time.values
ST = ds_silixa.st.values
AST = ds_silixa.ast.values

Now this data has to be inserted into an xarray Dataset

ds = xr.Dataset()
ds["x"] = ("x", x)
ds["time"] = ("time", time)
ds["st"] = (["x", "time"], ST)
ds["ast"] = (["x", "time"], AST)

ds = DataStore(ds)
print(ds)

<dtscalibration.DataStore>
Sections:                  ()
Dimensions:    (x: 1461, time: 3, trans_att: 0)
Coordinates:
  * x          (x) float64 -80.74 -80.62 -80.49 -80.36 ... 104.6 104.7 104.8
  * time       (time) datetime64[ns] 2018-05-04T12:22:17.710000 ... 2018-05-0...
  * trans_att  (trans_att) float64
Data variables:
    st         (x, time) float64 -0.8058 0.4287 -0.513 ... 27.99 27.83 28.81
    ast        (x, time) float64 -0.2459 -0.5932 0.1111 ... 36.2 35.7 35.16
Attributes:
    _sections:  null\n...\n

For calibration, a few more paramaters are needed:

• acquisition time (for calculating residuals for WLS calibration)

• reference temperatures

• a double ended flag

We’ll put these into the custom DataStore:

ds["acquisitiontimeFW"] = ds_silixa["acquisitiontimeFW"].values
ds["userAcquisitionTimeFW"] = ds_silixa["acquisitiontimeFW"].values
ds["temp1"] = ds_silixa["probe1Temperature"]
ds["temp2"] = ds_silixa["probe2Temperature"]

ds.attrs["isDoubleEnded"] = "0"

Now we can calibrate the data as usual (ordinary least squares in this example).

ds = ds.sel(x=slice(-30, 101))
sections = {
    "temp1": [slice(20, 25.5)],  # warm bath
    "temp2": [slice(5.5, 15.5)],  # cold bath
}
ds.sections = sections

st_var, resid = ds.variance_stokes_constant(st_label="st")
ast_var, _ = ds.variance_stokes_constant(st_label="ast")
ds.calibration_single_ended(st_var=st_var, ast_var=ast_var)

ds.isel(time=0).tmpf.plot();

ds

/home/docs/checkouts/readthedocs.org/user_builds/python-dts-calibration/envs/latest/lib/python3.9/site-packages/xarray/core/dataarray.py:861: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  return key in self.data
/home/docs/checkouts/readthedocs.org/user_builds/python-dts-calibration/envs/latest/lib/python3.9/site-packages/xarray/core/dataarray.py:861: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  return key in self.data

<dtscalibration.DataStore>
Sections:
    temp1                  ( 18.02 +/- 0.00°C)      20.00 - 25.50
    temp2                  (  6.62 +/- 0.00°C)      5.50 - 15.50
Dimensions:                (x: 1030, time: 3, trans_att: 0,
                            acquisitiontimeFW: 3, userAcquisitionTimeFW: 3,
                            comp_fw: 12, params1: 5, params2: 5)
Coordinates:
  * x                      (x) float64 -29.9 -29.78 -29.65 ... 100.6 100.8 100.9
  * time                   (time) datetime64[ns] 2018-05-04T12:22:17.710000 ....
  * trans_att              (trans_att) float64
  * acquisitiontimeFW      (acquisitiontimeFW) timedelta64[ns] 00:00:30 ... 0...
  * userAcquisitionTimeFW  (userAcquisitionTimeFW) timedelta64[ns] 00:00:30 ....
    filename               (time) <U31 'channel 2_20180504132202074.xml' ... ...
    filename_tstamp        (time) int64 20180504132202074 ... 20180504132303723
    timestart              (time) datetime64[ns] 2018-05-04T12:22:02.710000 ....
    timeend                (time) datetime64[ns] 2018-05-04T12:22:32.710000 ....
  * comp_fw                (comp_fw) object 'dT_dst' 'dT_dast' ... 'dta_ddalpha'
Dimensions without coordinates: params1, params2
Data variables: (12/21)
    st                     (x, time) float64 6.267e+03 6.272e+03 ... 2.619e+03
    ast                    (x, time) float64 5.473e+03 5.473e+03 ... 2.09e+03
    temp1                  (time) float32 18.02 18.02 18.02
    temp2                  (time) float32 6.62 6.617 6.617
    gamma                  float64 481.9
    gamma_var              float64 0.4209
    ...                     ...
    talpha_fw_full_var     (x, time) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    tmpf                   (x, time) float64 25.46 25.36 25.47 ... 10.1 10.18
    var_fw_da              (comp_fw, x, time) float64 0.0003519 ... 0.0
    tmpf_var               (x, time) float64 0.006567 0.006543 ... 0.01329
    p_val                  (params1) float64 481.9 -2.073e-05 1.478 1.477 1.477
    p_cov                  (params1, params2) float64 0.4209 ... 5.666e-06
Attributes:
    _sections:      temp1:\n- !!python/object/apply:builtins.slice\n  - 20.0\...
    isDoubleEnded:  0

xarray.DataStore

• Dimensions:
  • x: 1030
  • time: 3
  • trans_att: 0
  • acquisitiontimeFW: 3
  • userAcquisitionTimeFW: 3
  • comp_fw: 12
  • params1: 5
  • params2: 5
• Coordinates: (10)
  • x
    (x)
    float64
    -29.9 -29.78 -29.65 ... 100.8 100.9

    name :

    distance

    description :

    Length along fiber

    long_description :

    Starting at connector of forward channel

    units :

    m

    array([-29.9044, -29.7773, -29.6502, ..., 100.627 , 100.754 , 100.881 ])

  • time
    (time)
    datetime64[ns]
    2018-05-04T12:22:17.710000 ... 2...

    description :

    time halfway the measurement

    timezone :

    UTC

    array(['2018-05-04T12:22:17.710000000', '2018-05-04T12:22:47.702000000',
           '2018-05-04T12:23:18.716000000'], dtype='datetime64[ns]')

  • trans_att
    (trans_att)
    float64

    name :

    Locations introducing transient directional differential attenuation

    description :

    Locations along the x-dimension introducing transient directional differential attenuation

    long_description :

    Connectors introduce additional differential attenuation that is different for the forward and backward direction, and varies over time.

    units :

    m

    array([], dtype=float64)

  • acquisitiontimeFW
    (acquisitiontimeFW)
    timedelta64[ns]
    00:00:30 00:00:30 00:00:30

    array([30000000000, 30000000000, 30000000000], dtype='timedelta64[ns]')

  • userAcquisitionTimeFW
    (userAcquisitionTimeFW)
    timedelta64[ns]
    00:00:30 00:00:30 00:00:30

    array([30000000000, 30000000000, 30000000000], dtype='timedelta64[ns]')

  • filename
    (time)
    <U31
    'channel 2_20180504132202074.xml...

    array(['channel 2_20180504132202074.xml',
           'channel 2_20180504132232903.xml',
           'channel 2_20180504132303723.xml'], dtype='<U31')

  • filename_tstamp
    (time)
    int64
    20180504132202074 ... 2018050413...

    array([20180504132202074, 20180504132232903, 20180504132303723])

  • timestart
    (time)
    datetime64[ns]
    2018-05-04T12:22:02.710000 ... 2...

    description :

    time start of the measurement

    timezone :

    UTC

    array(['2018-05-04T12:22:02.710000000', '2018-05-04T12:22:32.702000000',
           '2018-05-04T12:23:03.716000000'], dtype='datetime64[ns]')

  • timeend
    (time)
    datetime64[ns]
    2018-05-04T12:22:32.710000 ... 2...

    description :

    time end of the measurement

    timezone :

    UTC

    array(['2018-05-04T12:22:32.710000000', '2018-05-04T12:23:02.702000000',
           '2018-05-04T12:23:33.716000000'], dtype='datetime64[ns]')

  • comp_fw
    (comp_fw)
    object
    'dT_dst' ... 'dta_ddalpha'

    array(['dT_dst', 'dT_dast', 'dT_gamma', 'dT_dc', 'dT_ddalpha', 'dT_dta',
           'dgamma_dc', 'dta_dgamma', 'dta_dc', 'dgamma_ddalpha', 'ddalpha_dc',
           'dta_ddalpha'], dtype=object)

• Data variables: (21)
  • st
    (x, time)
    float64
    6.267e+03 6.272e+03 ... 2.619e+03

    array([[6267.02, 6272.41, 6276.24],
           [6206.54, 6211.72, 6214.93],
           [6057.48, 6065.36, 6066.5 ],
           ...,
           [2634.95, 2623.17, 2618.16],
           [2634.59, 2624.87, 2617.84],
           [2635.04, 2624.33, 2618.71]])

  • ast
    (x, time)
    float64
    5.473e+03 5.473e+03 ... 2.09e+03

    array([[5473.15, 5473.18, 5479.61],
           [5386.9 , 5387.21, 5392.27],
           [5242.59, 5243.28, 5248.54],
           ...,
           [2102.93, 2094.09, 2090.71],
           [2102.88, 2093.12, 2090.87],
           [2103.14, 2093.44, 2089.98]])

  • temp1
    (time)
    float32
    18.02 18.02 18.02

    array([18.0204, 18.0211, 18.0216], dtype=float32)

  • temp2
    (time)
    float32
    6.62 6.617 6.617

    array([6.61986, 6.61692, 6.61695], dtype=float32)

  • gamma
    ()
    float64
    481.9

    array(481.91651701)

  • gamma_var
    ()
    float64
    0.4209

    array(0.42086891)

  • talpha_fw
    (trans_att, time)
    float64

    array([], shape=(0, 3), dtype=float64)

  • talpha_fw_var
    (trans_att, time)
    float64

    array([], shape=(0, 3), dtype=float64)

  • c
    (time)
    float64
    1.478 1.477 1.477

    array([1.47781366, 1.47746845, 1.47744336])

  • c_var
    (time)
    float64
    5.666e-06 5.666e-06 5.666e-06

    array([5.66556368e-06, 5.66580611e-06, 5.66557471e-06])

  • dalpha
    ()
    float64
    -2.073e-05

    array(-2.07252964e-05)

  • dalpha_var
    ()
    float64
    4.66e-11

    array(4.66027148e-11)

  • alpha
    (x)
    float64
    0.0006198 0.0006171 ... -0.002091

    array([ 0.00061978,  0.00061714,  0.00061451, ..., -0.00208552,
           -0.00208816, -0.00209079])

  • alpha_var
    (x)
    float64
    4.168e-08 4.132e-08 ... 4.743e-07

    array([4.16755561e-08, 4.13220491e-08, 4.09700479e-08, ...,
           4.71889450e-07, 4.73081332e-07, 4.74274718e-07])

  • talpha_fw_full
    (x, time)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0.],
           [0., 0., 0.],
           [0., 0., 0.],
           ...,
           [0., 0., 0.],
           [0., 0., 0.],
           [0., 0., 0.]])

  • talpha_fw_full_var
    (x, time)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0.],
           [0., 0., 0.],
           [0., 0., 0.],
           ...,
           [0., 0., 0.],
           [0., 0., 0.],
           [0., 0., 0.]])

  • tmpf
    (x, time)
    float64
    25.46 25.36 25.47 ... 10.1 10.18

    name :

    tmpf

    description :

    Temperature estimated using the Stokes and anti-Stokes from the Forward channel

    units :

    $^\circ$C

    array([[25.45738942, 25.36324094, 25.47215385],
           [24.3175255 , 24.23831562, 24.32042176],
           [23.79658962, 23.6460403 , 23.79963163],
           ...,
           [10.12022944, 10.22239788, 10.27596997],
           [10.13946049, 10.03780565, 10.3095433 ],
           [10.13204591, 10.09793338, 10.18365264]])

  • var_fw_da
    (comp_fw, x, time)
    float64
    0.0003519 0.0003508 ... 0.0 0.0

    array([[[ 0.00035188,  0.00035083,  0.00035092],
            [ 0.00035333,  0.00035236,  0.00035239],
            [ 0.00036834,  0.00036664,  0.00036726],
            ...,
            [ 0.00161204,  0.00162889,  0.00163637],
            [ 0.00161291,  0.00162255,  0.00163755],
            [ 0.00161219,  0.0016246 ,  0.00163355]],
    
           [[ 0.00035394,  0.00035349,  0.00035318],
            [ 0.00035982,  0.00035939,  0.00035911],
            [ 0.00037724,  0.00037638,  0.00037641],
            ...,
            [ 0.00194158,  0.00196083,  0.00196867],
            [ 0.0019422 ,  0.00195754,  0.0019693 ],
            [ 0.00194152,  0.00195861,  0.00196748]],
    
           [[ 0.16158622,  0.16148434,  0.1616022 ],
            [ 0.16035494,  0.16026955,  0.16035806],
            [ 0.15979379,  0.15963181,  0.15979707],
            ...,
    ...
            ...,
            [ 0.08011869,  0.08020541,  0.08025091],
            [ 0.08023614,  0.0801498 ,  0.08038075],
            [ 0.08033097,  0.08030196,  0.08037488]],
    
           [[ 0.03077049,  0.03073327,  0.03077693],
            [ 0.03017454,  0.03014395,  0.03017606],
            [ 0.02983583,  0.02977689,  0.02983739],
            ...,
            [-0.08385237, -0.0839777 , -0.08403789],
            [-0.083981  , -0.08386481, -0.08418383],
            [-0.08407805, -0.08404186, -0.0841403 ]],
    
           [[-0.        , -0.        , -0.        ],
            [-0.        , -0.        , -0.        ],
            [-0.        , -0.        , -0.        ],
            ...,
            [ 0.        ,  0.        ,  0.        ],
            [ 0.        ,  0.        ,  0.        ],
            [ 0.        ,  0.        ,  0.        ]]])

  • tmpf_var
    (x, time)
    float64
    0.006567 0.006543 ... 0.01329

    name :

    tmpf_var

    description :

    Uncertainty variance in tmpf estimated with linear-error propagation

    units :

    $^\circ$C$^2$

    array([[0.00656749, 0.00654302, 0.00656801],
           [0.0062431 , 0.00622356, 0.00624042],
           [0.00611552, 0.00607668, 0.0061126 ],
           ...,
           [0.01319432, 0.01321602, 0.01322631],
           [0.01322211, 0.01325438, 0.01325288],
           [0.01324987, 0.01327945, 0.01328885]])

  • p_val
    (params1)
    float64
    481.9 -2.073e-05 1.478 1.477 1.477

    array([ 4.81916517e+02, -2.07252964e-05,  1.47781366e+00,  1.47746845e+00,
            1.47744336e+00])

  • p_cov
    (params1, params2)
    float64
    0.4209 -4.068e-06 ... 5.666e-06

    array([[ 4.20868910e-01, -4.06752397e-06,  1.54382836e-03,
             1.54385882e-03,  1.54383948e-03],
           [-4.06752397e-06,  4.66027148e-11, -1.50283292e-08,
            -1.50290980e-08, -1.50285015e-08],
           [ 1.54382836e-03, -1.50283292e-08,  5.66556368e-06,
             5.66477469e-06,  5.66469771e-06],
           [ 1.54385882e-03, -1.50290980e-08,  5.66477469e-06,
             5.66580611e-06,  5.66481644e-06],
           [ 1.54383948e-03, -1.50285015e-08,  5.66469771e-06,
             5.66481644e-06,  5.66557471e-06]])

• Indexes: (6)
  • x
    PandasIndex

    PandasIndex(Index([-29.9044, -29.7773, -29.6502, -29.5231,  -29.396, -29.2689, -29.1418,
           -29.0147, -28.8876, -28.7605,
           ...
            99.7372,  99.8643,  99.9914,  100.119,  100.246,  100.373,    100.5,
            100.627,  100.754,  100.881],
          dtype='float64', name='x', length=1030))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2018-05-04 12:22:17.710000', '2018-05-04 12:22:47.702000',
                   '2018-05-04 12:23:18.716000'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • trans_att
    PandasIndex

    PandasIndex(Index([], dtype='float64', name='trans_att'))

  • acquisitiontimeFW
    PandasIndex

    PandasIndex(TimedeltaIndex(['0 days 00:00:30', '0 days 00:00:30', '0 days 00:00:30'], dtype='timedelta64[ns]', name='acquisitiontimeFW', freq=None))

  • userAcquisitionTimeFW
    PandasIndex

    PandasIndex(TimedeltaIndex(['0 days 00:00:30', '0 days 00:00:30', '0 days 00:00:30'], dtype='timedelta64[ns]', name='userAcquisitionTimeFW', freq=None))

  • comp_fw
    PandasIndex

    PandasIndex(Index(['dT_dst', 'dT_dast', 'dT_gamma', 'dT_dc', 'dT_ddalpha', 'dT_dta',
           'dgamma_dc', 'dta_dgamma', 'dta_dc', 'dgamma_ddalpha', 'ddalpha_dc',
           'dta_ddalpha'],
          dtype='object', name='comp_fw'))

• Attributes: (2)

  _sections :

  temp1: - !!python/object/apply:builtins.slice - 20.0 - 25.5 - null temp2: - !!python/object/apply:builtins.slice - 5.5 - 15.5 - null

  isDoubleEnded :

  0