search

sappelhoff / pyprep

A Python implementation of the Preprocessing Pipeline (PREP) for EEG data
https://pyprep.readthedocs.io/en/latest/
MIT License
128 stars 30 forks source link

Saving prep.raw give error " ValueError: Measurement infos are inconsistent for dig" #113

Closed walkerped closed 2 years ago

walkerped commented 2 years ago

I am trying to data from openneuro (the first subject under https://openneuro.org/datasets/ds003690/versions/1.0.0) through pyprep with the following code.

import os
import mne
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
from pyprep.prep_pipeline import PrepPipeline

subs = ['10']
data_dir = '/mnt/f/portfolio/data/restEEG'
study_dir = os.path.join(data_dir,'ds003690')

for sub in subs:

    # load file
    subStr = f'sub-AB{sub}'
    sub_dir = os.path.join(study_dir,subStr,'eeg')
    eeg_file = os.path.join(sub_dir,f'sub-AB{sub}_task-passive_run-1_eeg.set')
    raw_eeg = mne.io.read_raw_eeglab(eeg_file)

    # Set ecg and eog
    raw_eeg = raw_eeg.set_channel_types({'EKG':'ecg','VEO':'eog','HEO':'eog' })
    raw_eeg.drop_channels(['CB1','CB2','R-Dia-X-(mm)','R-Dia-Y-(mm)'])

    # Extract sample rate
    sample_rate = raw_eeg.info["sfreq"]

    # Make a copy of the data
    prep_eeg = raw_eeg.copy()

    # Run prep
    prep_params = {
    "ref_chs": "eeg",
    "reref_chs": "eeg",
    "line_freqs": np.arange(60, sample_rate / 2, 60),
    }
    prep = PrepPipeline(prep_eeg, prep_params, mont1020)
    prep.fit()

The command appears to be successful giving the following output:

/tmp/ipykernel_16702/2558830226.py:9: RuntimeWarning: Unknown types found, setting as type EEG:
heog: ['HEO']
pupil: ['R-Dia-X-(mm)', 'R-Dia-Y-(mm)']
veog: ['VEO']
  raw_eeg = mne.io.read_raw_eeglab(eeg_file)
/tmp/ipykernel_16702/2558830226.py:9: RuntimeWarning: Data will be preloaded. preload=False or a string preload is not supported when the data is stored in the .set file
  raw_eeg = mne.io.read_raw_eeglab(eeg_file)
/tmp/ipykernel_16702/2558830226.py:9: RuntimeWarning: Not setting position of 1 ecg channel found in montage:
['EKG']
Consider setting the channel types to be of EEG/sEEG/ECoG/DBS/fNIRS using inst.set_channel_types before calling inst.set_montage, or omit these channels when creating your montage.
  raw_eeg = mne.io.read_raw_eeglab(eeg_file)

Setting up high-pass filter at 1 Hz

FIR filter parameters

Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1651 samples (3.302 sec)

Setting up high-pass filter at 1 Hz

FIR filter parameters

Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1651 samples (3.302 sec)

Removed notch frequencies (Hz):
     60.00 : 3127 windows
    120.00 : 3127 windows
    180.00 : 3127 windows
    239.00 :   59 windows
    240.00 : 3127 windows
    241.00 :   59 windows
Setting up high-pass filter at 1 Hz

FIR filter parameters

Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1651 samples (3.302 sec)

Executing RANSAC
This may take a while, so be patient...
Progress: 10%... 20%... 30%... 40%... 50%... 60%... 70%... 80%... 90%... 

2022-03-26 15:16:31,717 - pyprep.reference - INFO - Bad channels: {'bad_by_nan': [], 'bad_by_flat': [], 'bad_by_deviation': [], 'bad_by_hf_noise': ['Cz', 'CPz'], 'bad_by_correlation': ['Cz', 'CPz'], 'bad_by_SNR': ['Cz', 'CPz'], 'bad_by_dropout': [], 'bad_by_ransac': [], 'bad_all': ['Cz', 'CPz']}

100%

RANSAC done!
Executing RANSAC
This may take a while, so be patient...
Progress: 10%... 20%... 30%... 40%... 50%... 60%... 70%... 80%... 90%... 

2022-03-26 15:16:52,691 - pyprep.reference - INFO - Bad channels: {'bad_by_nan': [], 'bad_by_flat': [], 'bad_by_deviation': [], 'bad_by_hf_noise': ['Cz', 'CPz'], 'bad_by_correlation': ['AF4', 'Cz', 'CPz'], 'bad_by_SNR': ['Cz', 'CPz'], 'bad_by_dropout': [], 'bad_by_ransac': ['M1'], 'bad_all': ['AF4', 'Cz', 'M1', 'CPz']}

100%

RANSAC done!
Interpolating bad channels
    Automatic origin fit: head of radius 95.5 mm
Computing interpolation matrix from 55 sensor positions
Interpolating 4 sensors

2022-03-26 15:16:52,939 - pyprep.reference - INFO - Iterations: 1

Executing RANSAC
This may take a while, so be patient...
Progress: 10%... 20%... 30%... 40%... 50%... 60%... 70%... 80%... 90%... 

2022-03-26 15:17:07,594 - pyprep.reference - INFO - Bad channels: {'bad_by_nan': [], 'bad_by_flat': [], 'bad_by_deviation': [], 'bad_by_hf_noise': ['Cz', 'CPz'], 'bad_by_correlation': ['AF4', 'Cz', 'CPz'], 'bad_by_SNR': ['Cz', 'CPz'], 'bad_by_dropout': [], 'bad_by_ransac': ['M1'], 'bad_all': ['AF4', 'Cz', 'M1', 'CPz']}

100%

RANSAC done!
Interpolating bad channels
    Automatic origin fit: head of radius 95.5 mm
Computing interpolation matrix from 55 sensor positions
Interpolating 4 sensors

2022-03-26 15:17:07,844 - pyprep.reference - INFO - Iterations: 2

Executing RANSAC
This may take a while, so be patient...
Progress: 10%... 20%... 30%... 40%... 50%... 60%... 70%... 80%... 90%... 

2022-03-26 15:17:22,721 - pyprep.reference - INFO - Bad channels: {'bad_by_nan': [], 'bad_by_flat': [], 'bad_by_deviation': [], 'bad_by_hf_noise': ['Cz', 'CPz'], 'bad_by_correlation': ['AF4', 'Cz', 'CPz'], 'bad_by_SNR': ['Cz', 'CPz'], 'bad_by_dropout': [], 'bad_by_ransac': ['M1'], 'bad_all': ['AF4', 'Cz', 'M1', 'CPz']}
2022-03-26 15:17:22,722 - pyprep.reference - INFO - Robust reference done

100%

RANSAC done!
Interpolating bad channels
    Automatic origin fit: head of radius 95.5 mm
Computing interpolation matrix from 55 sensor positions
Interpolating 4 sensors
Setting up high-pass filter at 1 Hz

FIR filter parameters

Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1651 samples (3.302 sec)

Executing RANSAC
This may take a while, so be patient...
Progress: 10%... 20%... 30%... 40%... 50%... 60%... 70%... 80%... 90%... 100%

RANSAC done!
Found 1 uniquely bad channels:

0 by NaN: []

0 by flat: []

0 by deviation: []

0 by HF noise: []

1 by correlation: ['AF4']

0 by SNR: []

0 by dropout: []

0 by RANSAC: []

Interpolating bad channels
    Automatic origin fit: head of radius 95.5 mm
Computing interpolation matrix from 56 sensor positions
Interpolating 3 sensors
Setting up high-pass filter at 1 Hz

FIR filter parameters

Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1651 samples (3.302 sec)

Executing RANSAC
This may take a while, so be patient...
Progress: 10%... 20%... 30%... 40%... 50%... 60%... 70%... 80%... 90%... 100%

RANSAC done!

And I am able to retrieve info about bad channels with the following code and output:

print("Bad channels original: {}".format(prep.noisy_channels_original["bad_all"]))
print("Bad channels after interpolation: {}".format(prep.still_noisy_channels))

Bad channels: ['AF4', 'Cz', 'CPz']
Bad channels original: ['Cz', 'CPz']
Bad channels after interpolation: []

But when I try to refer to prep.raw in any command, for example EEG = prep.raw.get_data(), I get an error like:

ValueError                                Traceback (most recent call last)
Input In [106], in <cell line: 1>()
----> 1 EEG = prep.raw.get_data()

File ~/anaconda3/envs/restEEG/lib/python3.9/site-packages/pyprep/prep_pipeline.py:175, in PrepPipeline.raw(self)
    173     return full_raw
    174 else:
--> 175     return full_raw.add_channels([self.raw_non_eeg])

File ~/anaconda3/envs/restEEG/lib/python3.9/site-packages/mne/channels/channels.py:954, in UpdateChannelsMixin.add_channels(self, add_list, force_update_info)
    952 # Create final data / info objects
    953 infos = [self.info] + [inst.info for inst in add_list]
--> 954 new_info = _merge_info(infos, force_update_to_first=force_update_info)
    956 # Now update the attributes
    957 if isinstance(self._data, np.memmap) and con_axis == 0 and \
    958         sys.platform != 'darwin':  # resizing not available--no mremap
    959     # Use a resize and fill in other ones

File <decorator-gen-35>:12, in _merge_info(infos, force_update_to_first, verbose)

File ~/anaconda3/envs/restEEG/lib/python3.9/site-packages/mne/io/meas_info.py:2404, in _merge_info(infos, force_update_to_first, verbose)
   2402     else:
   2403         msg = ("Measurement infos are inconsistent for %s" % k)
-> 2404         raise ValueError(msg)
   2406 # other fields
   2407 other_fields = ['acq_pars', 'acq_stim', 'bads',
   2408                 'comps', 'custom_ref_applied', 'description',
   2409                 'experimenter', 'file_id', 'highpass', 'utc_offset',
   (...)
   2412                 'proj_id', 'proj_name', 'projs', 'sfreq', 'gantry_angle',
   2413                 'subject_info', 'sfreq', 'xplotter_layout', 'proc_history']

ValueError: Measurement infos are inconsistent for dig

Although other attributes of prep seem to be fine:

prep.EEG_new

array([[ 6.92872951e-19, -2.80886176e-06, -5.02157417e-06, ...,
         5.30101354e-06,  2.26315846e-06,  5.75982404e-19],
       [ 1.00119294e-18, -2.10187357e-06, -2.50465512e-06, ...,
         2.35277111e-07, -3.63729870e-06, -7.86046575e-19],
       [-7.64950161e-19,  1.81967321e-06,  1.10303559e-05, ...,
         1.34994738e-05,  7.49501261e-06, -6.28498447e-19],
       ...,
       [ 4.10578045e-18, -2.07593091e-07,  6.23457305e-08, ...,
         8.15795696e-06,  1.91804891e-06,  1.62630326e-19],
       [-1.05879118e-18, -6.16081013e-07, -7.84905784e-07, ...,
         8.63973602e-06,  1.81714657e-06,  5.14996032e-18],
       [-3.32936888e-18, -1.16585488e-06, -1.58660128e-06, ...,
         9.17628802e-06,  3.79875150e-06,  1.15196481e-18]])
sappelhoff commented 2 years ago

Thanks for the report, this might have to do with the recent MNE-Python 1.0 release ... are you using mne at version 1.0?

Perhaps the warning is spurious for our case and could be suppressed.

walkerped commented 2 years ago

Yes, I am using MNE v. 1.0.0. It seems to be causing any process I run with prep.raw to fail, rather than just being a warning. So if I run EEG = prep.raw.get_data() and get the error, I'll get a NameError: name 'EEG' is not defined if I try to refer to EEG. Similarly if I try to do prep.raw.plot(), I get the ValueError: Measurement infos are inconsistent for dig error again.

sappelhoff commented 2 years ago

please see what I wrote in #114

Also, please double check your code and see the warnings that your code produces:

/tmp/ipykernel_16702/2558830226.py:9: RuntimeWarning: Unknown types found, setting as type EEG: heog: ['HEO'] pupil: ['R-Dia-X-(mm)', 'R-Dia-Y-(mm)'] veog: ['VEO'] raw_eeg = mne.io.read_raw_eeglab(eeg_file) /tmp/ipykernel_16702/2558830226.py:9: RuntimeWarning: Data will be preloaded. preload=False or a string preload is not supported when the data is stored in the .set file raw_eeg = mne.io.read_raw_eeglab(eeg_file) /tmp/ipykernel_16702/2558830226.py:9: RuntimeWarning: Not setting position of 1 ecg channel found in montage: ['EKG'] Consider setting the channel types to be of EEG/sEEG/ECoG/DBS/fNIRS using inst.set_channel_types before calling inst.set_montage, or omit these channels when creating your montage. raw_eeg = mne.io.read_raw_eeglab(eeg_file)

walkerped commented 2 years ago

The update in https://github.com/sappelhoff/pyprep/pull/114 seems to have fixed it - I can plot prep.raw, and see that the bad channel was interpolated. Thanks!