MRG: Add support for GSR and temperature channels (#1059)

* Add support for GSR and temperature channels Also updates tiny_bids Will fail until mne-tools/mne-python#11108 has been merged Fixes #1046 * Indent * Preserve info * Require MNE devel to run tests * Fix unit writing * Write correct BIDS temperature unit (oC) * Better code * Add conditional * No warning * Skip test on MNE stable * Decorated the wrong test * Add channel statuses * Update doc/whats_new.rst Co-authored-by: Stefan Appelhoff <[email protected]> * Fix channel description * Manually curate tiny_bids to keep the diff small * Revert changes related to _orig_units handling * Fix test Co-authored-by: Stefan Appelhoff <[email protected]>
mne-tools · Aug 31, 2022 · 388616f · 388616f
1 parent c2f5ca1
commit 388616f
Show file tree

Hide file tree

Showing 11 changed files with 128 additions and 48 deletions.
diff --git a/doc/whats_new.rst b/doc/whats_new.rst
@@ -53,6 +53,8 @@ Detailed list of changes
 
 - :class:`~mne_bids.BIDSPath` now supports the BIDS "description" entity ``desc``, used in derivative data, by `Richard Höchenberger`_ (:gh:`1049`)
 
+- Added support for ``GSR`` (galvanic skin response / electrodermal activity, EDA) and ``TEMP`` (temperature) channel types, by `Richard Höchenberger`_ (:gh:`xxx`)
+
 🧐 API and behavior changes
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/mne_bids/config.py b/mne_bids/config.py
@@ -31,14 +31,18 @@
     '.sqd': 'KitYokogawa',
 }
 
-UNITS = {
+EXT_TO_UNIT_MAP = {
     '.con': 'm',
     '.ds': 'cm',
     '.fif': 'm',
     '.pdf': 'm',
     '.sqd': 'm'
 }
 
+UNITS_MNE_TO_BIDS_MAP = {
+    'C': 'oC',  # temperature in deg. C
+}
+
 meg_manufacturers = {
     '.con': 'KIT/Yokogawa',
     '.ds': 'CTF',

diff --git a/mne_bids/tests/data/tiny_bids/code/make_tiny_bids_dataset.py b/mne_bids/tests/data/tiny_bids/code/make_tiny_bids_dataset.py
@@ -1,8 +1,7 @@
 """Code used to generate the tiny_bids dataset."""
+
 # %%
 import json
-import os
-import os.path as op
 from pathlib import Path
 
 import mne
@@ -12,17 +11,20 @@
 from mne_bids import BIDSPath, write_raw_bids
 
 data_path = mne.datasets.testing.data_path()
-vhdr_fname = op.join(data_path, "montage", "bv_dig_test.vhdr")
-captrak_path = op.join(data_path, "montage", "captrak_coords.bvct")
+vhdr_path = data_path / "montage" / "bv_dig_test.vhdr"
+captrak_path = data_path / "montage" / "captrak_coords.bvct"
 
-mne_bids_root = Path(mne_bids.__file__).parent.parent
-tiny_bids = op.join(mne_bids_root, "mne_bids", "tests", "data", "tiny_bids")
-os.makedirs(tiny_bids, exist_ok=True)
+mne_bids_root = Path(mne_bids.__file__).parents[1]
+tiny_bids_root = mne_bids_root / "mne_bids" / "tests" / "data" / "tiny_bids"
+tiny_bids_root.mkdir(exist_ok=True)
 
-bids_path = BIDSPath(subject="01", task="rest", session="eeg", root=tiny_bids)
+bids_path = BIDSPath(
+    subject="01", task="rest", session="eeg", suffix="eeg", extension=".vhdr",
+    datatype="eeg", root=tiny_bids_root
+)
 
 # %%
-raw = mne.io.read_raw_brainvision(vhdr_fname)
+raw = mne.io.read_raw_brainvision(vhdr_path, preload=True)
 montage = mne.channels.read_dig_captrak(captrak_path)
 
 raw.set_channel_types(dict(ECG="ecg", HEOG="eog", VEOG="eog"))
@@ -39,25 +41,67 @@
     "hand": 3,
 }
 
-raw.set_annotations(None)
-events = np.array([[0, 0, 1], [1000, 0, 2]])
-event_id = {"start_experiment": 1, "show_stimulus": 2}
+# %%
+# Add GSR and temperature channels
+if 'GSR' not in raw.ch_names and 'Temperature' not in raw.ch_names:
+    gsr_data = np.array([2.1e-6] * len(raw.times))
+    temperature_data = np.array([36.5] * len(raw.times))
+
+    gsr_and_temp_data = np.concatenate([
+        np.atleast_2d(gsr_data),
+        np.atleast_2d(temperature_data),
+    ])
+    gsr_and_temp_info = mne.create_info(
+        ch_names=["GSR", "Temperature"],
+        sfreq=raw.info["sfreq"],
+        ch_types=["gsr", "temperature"],
+    )
+    gsr_and_temp_info["line_freq"] = raw.info["line_freq"]
+    gsr_and_temp_info["subject_info"] = raw.info["subject_info"]
+    with gsr_and_temp_info._unlock():
+        gsr_and_temp_info["lowpass"] = raw.info["lowpass"]
+        gsr_and_temp_info["highpass"] = raw.info["highpass"]
+    gsr_and_temp_raw = mne.io.RawArray(
+        data=gsr_and_temp_data,
+        info=gsr_and_temp_info,
+        first_samp=raw.first_samp,
+    )
+    raw.add_channels([gsr_and_temp_raw])
+    del gsr_and_temp_raw, gsr_and_temp_data, gsr_and_temp_info
 
 # %%
+raw.set_annotations(None)
+events = np.array([
+    [0, 0, 1],
+    [1000, 0, 2]
+])
+event_id = {
+    "start_experiment": 1,
+    "show_stimulus": 2
+}
 
+# %%
 write_raw_bids(
-    raw, bids_path, events=events, event_id=event_id, overwrite=True
+    raw, bids_path, events=events, event_id=event_id, overwrite=True,
+    allow_preload=True, format="BrainVision",
+)
+mne_bids.mark_channels(
+    bids_path=bids_path,
+    ch_names=['C3', 'C4', 'PO10', 'GSR', 'Temperature'],
+    status=['good', 'good', 'bad', 'good', 'good'],
+    descriptions=['resected', 'resected', 'continuously flat',
+                  'left index finger', 'left ear']
 )
 
 # %%
-
-dataset_description_json = op.join(tiny_bids, "dataset_description.json")
-with open(dataset_description_json, "r") as fin:
-    ds_json = json.load(fin)
+dataset_description_json_path = tiny_bids_root / "dataset_description.json"
+ds_json = json.loads(
+    dataset_description_json_path.read_text(encoding="utf-8")
+)
 
 ds_json["Name"] = "tiny_bids"
 ds_json["Authors"] = ["MNE-BIDS Developers", "And Friends"]
 
-with open(dataset_description_json, "w") as fout:
+with open(dataset_description_json_path, "w", encoding="utf-8") as fout:
     json.dump(ds_json, fout, indent=4)
     fout.write("\n")
diff --git a/mne_bids/tests/data/tiny_bids/sub-01/ses-eeg/eeg/sub-01_ses-eeg_electrodes.tsv b/mne_bids/tests/data/tiny_bids/sub-01/ses-eeg/eeg/sub-01_ses-eeg_electrodes.tsv
@@ -63,6 +63,3 @@ PO3	-0.03295276056603997	-0.09836112961335748	0.09818535662931108	0.0
 POz	0.00141250247182117	-0.09659898488597438	0.1092860795662283	2.0
 PO4	0.03732555069518968	-0.09330847523481997	0.09871963971555255	0.0
 PO8	0.061422359682971855	-0.08824937021900595	0.06730156974139692	0.0
-ECG	n/a	n/a	n/a	n/a
-HEOG	n/a	n/a	n/a	n/a
-VEOG	n/a	n/a	n/a	n/a
diff --git a/mne_bids/tests/data/tiny_bids/sub-01/ses-eeg/eeg/sub-01_ses-eeg_task-rest_channels.tsv b/mne_bids/tests/data/tiny_bids/sub-01/ses-eeg/eeg/sub-01_ses-eeg_task-rest_channels.tsv
@@ -66,3 +66,5 @@ PO8	EEG	µV	0.015915494309189534	1000.0	ElectroEncephaloGram	5000.0	good	n/a
 ECG	ECG	µV	0.015915494309189534	1000.0	ElectroCardioGram	5000.0	good	n/a
 HEOG	EOG	µV	0.015915494309189534	1000.0	ElectroOculoGram	5000.0	good	n/a
 VEOG	EOG	µV	0.015915494309189534	1000.0	ElectroOculoGram	5000.0	good	n/a
+GSR	GSR	S	0.015915494309189534	1000.0	Galvanic skin response (electrodermal activity, EDA)	5000.0	good	left index finger
+Temperature	TEMP	oC	0.015915494309189534	1000.0	Temperature	5000.0	good	left ear
diff --git a/mne_bids/tests/data/tiny_bids/sub-01/ses-eeg/eeg/sub-01_ses-eeg_task-rest_eeg.eeg b/mne_bids/tests/data/tiny_bids/sub-01/ses-eeg/eeg/sub-01_ses-eeg_task-rest_eeg.eeg
diff --git a/mne_bids/tests/data/tiny_bids/sub-01/ses-eeg/eeg/sub-01_ses-eeg_task-rest_eeg.vhdr b/mne_bids/tests/data/tiny_bids/sub-01/ses-eeg/eeg/sub-01_ses-eeg_task-rest_eeg.vhdr
@@ -8,7 +8,7 @@ MarkerFile=sub-01_ses-eeg_task-rest_eeg.vmrk
 DataFormat=BINARY
 ; Data orientation: MULTIPLEXED=ch1,pt1, ch2,pt1 ...
 DataOrientation=MULTIPLEXED
-NumberOfChannels=67
+NumberOfChannels=69
 ; Sampling interval in microseconds
 SamplingInterval=200
 
@@ -87,6 +87,8 @@ Ch64=PO8,,0.1,µV
 Ch65=ECG,,0.1,µV
 Ch66=HEOG,,0.1,µV
 Ch67=VEOG,,0.1,µV
+Ch68=GSR,,0.1,S
+Ch69=Temperature,,0.1,C
 
 [Comment]
 
@@ -95,7 +97,7 @@ BrainVision Recorder Professional   -   V. 1.21.0303
 
 A m p l i f i e r  S e t u p
 ============================
-Number of channels: 67
+Number of channels: 69
 Sampling Rate [Hz]: 1000
 Sampling Interval [µS]: 1000
 
@@ -169,6 +171,8 @@ Channels
 65    ECG         65               0.1 µV             10             1000              Off                0
 66    HEOG        66               0.1 µV             10             1000              Off                0
 67    VEOG        67               0.1 µV             10             1000              Off                0
+68    GSR         68               0.1 µV             10             1000              Off                0
+69    Temperature 69               0.1 µV             10             1000              Off                0
 
 S o f t w a r e  F i l t e r s
 ==============================

diff --git a/mne_bids/tests/test_path.py b/mne_bids/tests/test_path.py
@@ -145,22 +145,23 @@ def test_search_folder_for_text(capsys):
     test_dir = op.dirname(__file__)
     search_folder_for_text('n/a', test_dir)
     captured = capsys.readouterr()
-    assert 'sub-01_ses-eeg_electrodes.tsv' in captured.out
-    assert '    1    name      x         y         z         impedance' in \
-        captured.out.split('\n')
-    assert '    66   ECG       n/a       n/a       n/a       n/a' in \
-        captured.out.split('\n')
-
+    assert 'sub-01_ses-eeg_task-rest_eeg.json' in captured.out
+    assert (
+        '    1    name      type      units     low_cutof high_cuto descripti sampling_ status    status_de\n'  # noqa: E501
+        '    2    Fp1       EEG       µV        0.0159154 1000.0    ElectroEn 5000.0    good      n/a'          # noqa: E501
+    ) in captured.out
     # test if pathlib.Path object
     search_folder_for_text('n/a', Path(test_dir))
 
     # test returning a string and without line numbers
     out = search_folder_for_text(
-        'n/a', test_dir, line_numbers=False, return_str=True)
-    assert 'sub-01_ses-eeg_electrodes.tsv' in out
-    assert '    name      x         y         z         impedance' in \
-        out.split('\n')
-    assert '    ECG       n/a       n/a       n/a       n/a' in out.split('\n')
+        'n/a', test_dir, line_numbers=False, return_str=True
+    )
+    assert 'sub-01_ses-eeg_task-rest_eeg.json' in out
+    assert (
+        '    name      type      units     low_cutof high_cuto descripti sampling_ status    status_de\n'  # noqa: E501
+        '    Fp1       EEG       µV        0.0159154 1000.0    ElectroEn 5000.0    good      n/a'          # noqa: E501
+    ) in out
 
 
 def test_print_dir_tree(capsys):

diff --git a/mne_bids/tests/test_read.py b/mne_bids/tests/test_read.py
@@ -58,8 +58,8 @@
 raw_fname_chpi = op.join(data_path, 'SSS', 'test_move_anon_raw.fif')
 
 # Tiny BIDS testing dataset
-mne_bids_root = Path(mne_bids.__file__).parent.parent
-tiny_bids = op.join(mne_bids_root, "mne_bids", "tests", "data", "tiny_bids")
+mne_bids_root = Path(mne_bids.__file__).parents[1]
+tiny_bids_root = mne_bids_root / "mne_bids" / "tests" / "data" / "tiny_bids"
 
 warning_str = dict(
     channel_unit_changed='ignore:The unit for chann*.:RuntimeWarning:mne',
@@ -572,6 +572,8 @@ def test_handle_scans_reading(tmp_path):
     assert new_acq_time != raw_01.info['meas_date']
 
 
+@requires_version('mne', '1.2')  # tiny_bids contains GSR & temperature chans
+@pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
 def test_handle_scans_reading_brainvision(tmp_path):
     """Test stability of BrainVision's different file extensions"""
     test_scan_eeg = OrderedDict(
@@ -591,9 +593,9 @@ def test_handle_scans_reading_brainvision(tmp_path):
         _to_tsv(test_scan, tmp_path / test_scan['filename'][0])
 
     bids_path = BIDSPath(subject='01', session='eeg', task='rest',
-                         datatype='eeg', root=tiny_bids)
-    with pytest.warns(RuntimeWarning, match='Not setting positions'):
-        raw = read_raw_bids(bids_path)
+                         datatype='eeg', root=tiny_bids_root)
+
+    raw = read_raw_bids(bids_path)
 
     for test_scan in [test_scan_eeg, test_scan_vmrk]:
         _handle_scans_reading(tmp_path / test_scan['filename'][0],
@@ -690,7 +692,6 @@ def test_handle_info_reading(tmp_path):
         assert raw.info['line_freq'] == 55
 
 
-@requires_version('mne', '0.24')
 @pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
 @pytest.mark.filterwarnings(warning_str['maxshield'])
 def test_handle_chpi_reading(tmp_path):
@@ -1272,3 +1273,16 @@ def test_file_not_found(tmp_path):
     bp.extension = None
     with pytest.raises(FileNotFoundError, match='File does not exist'):
         read_raw_bids(bids_path=bp)
+
+
+@requires_version('mne', '1.2')
+@pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
+def test_gsr_and_temp_reading():
+    """Test GSR and temperature channels are handled correctly."""
+    bids_path = BIDSPath(
+        subject='01', session='eeg', task='rest', datatype='eeg',
+        root=tiny_bids_root
+    )
+    raw = read_raw_bids(bids_path)
+    assert raw.get_channel_types(['GSR']) == ['gsr']
+    assert raw.get_channel_types(['Temperature']) == ['temperature']
diff --git a/mne_bids/utils.py b/mne_bids/utils.py
@@ -66,7 +66,8 @@ def _get_ch_type_mapping(fro='mne', to='bids'):
     if fro == 'mne' and to == 'bids':
         mapping = dict(eeg='EEG', misc='MISC', stim='TRIG', emg='EMG',
                        ecog='ECOG', seeg='SEEG', eog='EOG', ecg='ECG',
-                       resp='RESP', bio='MISC', dbs='DBS',
+                       resp='RESP', bio='MISC', dbs='DBS', gsr='GSR',
+                       temperature='TEMP',
                        # NIRS
                        fnirs_cw_amplitude='NIRSCWAMPLITUDE',
                        # MEG channels
@@ -78,7 +79,7 @@ def _get_ch_type_mapping(fro='mne', to='bids'):
     elif fro == 'bids' and to == 'mne':
         mapping = dict(EEG='eeg', MISC='misc', TRIG='stim', EMG='emg',
                        ECOG='ecog', SEEG='seeg', EOG='eog', ECG='ecg',
-                       RESP='resp',
+                       RESP='resp', GSR='gsr', TEMP='temperature',
                        # NIRS
                        NIRSCWAMPLITUDE='fnirs_cw_amplitude',
                        NIRS='fnirs_cw_amplitude',

diff --git a/mne_bids/write.py b/mne_bids/write.py
@@ -51,12 +51,13 @@
 from mne_bids.read import _find_matching_sidecar, _read_events
 from mne_bids.sidecar_updates import update_sidecar_json
 
-from mne_bids.config import (ORIENTATION, UNITS, MANUFACTURERS,
+from mne_bids.config import (ORIENTATION, EXT_TO_UNIT_MAP, MANUFACTURERS,
                              IGNORED_CHANNELS, ALLOWED_DATATYPE_EXTENSIONS,
                              BIDS_VERSION, REFERENCES, _map_options, reader,
                              ALLOWED_INPUT_EXTENSIONS, CONVERT_FORMATS,
                              ANONYMIZED_JSON_KEY_WHITELIST, PYBV_VERSION,
-                             BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS)
+                             BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS,
+                             UNITS_MNE_TO_BIDS_MAP,)
 
 
 _FIFF_SPLIT_SIZE = '2GB'  # MNE-Python default; can be altered during debugging
@@ -102,7 +103,10 @@ def _channels_tsv(raw, fname, overwrite=False):
                     ias='Internal Active Shielding',
                     dbs='Deep Brain Stimulation',
                     fnirs_cw_amplitude='Near Infrared Spectroscopy '
-                                       '(continuous wave)',)
+                                       '(continuous wave)',
+                    resp='Respiration',
+                    gsr='Galvanic skin response (electrodermal activity, EDA)',
+                    temperature='Temperature',)
     get_specific = ('mag', 'ref_meg', 'grad')
 
     # get the manufacturer from the file in the Raw object
@@ -120,12 +124,19 @@ def _channels_tsv(raw, fname, overwrite=False):
         ch_type.append(map_chs[_channel_type])
         description.append(map_desc[_channel_type])
     low_cutoff, high_cutoff = (raw.info['highpass'], raw.info['lowpass'])
-    if raw._orig_units:
+    if raw._orig_units is not None:
         units = [raw._orig_units.get(ch, 'n/a') for ch in raw.ch_names]
     else:
         units = [_unit2human.get(ch_i['unit'], 'n/a')
                  for ch_i in raw.info['chs']]
         units = [u if u not in ['NA'] else 'n/a' for u in units]
+
+    # Translate from MNE to BIDS unit naming
+    for idx, mne_unit in enumerate(units):
+        if mne_unit in UNITS_MNE_TO_BIDS_MAP:
+            bids_unit = UNITS_MNE_TO_BIDS_MAP[mne_unit]
+            units[idx] = bids_unit
+
     n_channels = raw.info['nchan']
     sfreq = raw.info['sfreq']
 
@@ -1764,7 +1775,7 @@ def write_raw_bids(
                 bids_path.update(extension='.vhdr')
     # Read in Raw object and extract metadata from Raw object if needed
     orient = ORIENTATION.get(ext, 'n/a')
-    unit = UNITS.get(ext, 'n/a')
+    unit = EXT_TO_UNIT_MAP.get(ext, 'n/a')
     manufacturer = MANUFACTURERS.get(ext, 'n/a')
 
     # save readme file unless it already exists