Exclude GDF channels (#8520)

* Exclude GDF channels

* Fix units

* Add test

* Add changelog entry

* Fix PEP8

* Fix GDF2 units

doc/changes/0.21.inc

@@ -56,6 +56,8 @@ Bugs
 
 - When reading BrainVision raw data, the channel units and types were sometimes not inferred correctly by `Richard Höchenberger`_ (:gh:`8434`)
 
+- Fix reading GDF files with excluded channels in :func:`mne.io.read_raw_gdf` (:gh:`8520` by `Clemens Brunner`_)
+
 .. _changes_0_21:
 
 Version 0.21
@@ -263,7 +265,7 @@ Enhancements
 
 - Add reader for SNIRF NIRS data in :func:`mne.io.read_raw_snirf` by `Robert Luke`_
 
-- `~mne.Evoked` has gained ``tmin`` and ``tmax`` attributes for more consistency with `~mne.Epochs` by `Richard Höchenberger`_		
+- `~mne.Evoked` has gained ``tmin`` and ``tmax`` attributes for more consistency with `~mne.Epochs` by `Richard Höchenberger`_
 
 Bugs
 ~~~~
@@ -448,7 +450,7 @@ API changes
 
   1. Arguments to clustering functions, such as `mne.stats.permutation_cluster_test`, and
   2. Function names for defining adjacency, such as `mne.spatio_temporal_src_adjacency` replacing ``mne.spatio_temporal_src_connectivity``.
-  
+
   The complete list of changed function names is:
 
   - ``mne.channels.find_ch_connectivity`` → `~mne.channels.find_ch_adjacency`

mne/io/edf/edf.py

@@ -719,18 +719,20 @@ def _read_gdf_header(fname, exclude):
             channels = list(range(nchan))
             ch_names = [fid.read(16).decode('latin-1').strip(' \x00')
                         for ch in channels]
+            exclude = _find_exclude_idx(ch_names, exclude)
+            sel = np.setdiff1d(np.arange(len(ch_names)), exclude)
             fid.seek(80 * len(channels), 1)  # transducer
             units = [fid.read(8).decode('latin-1').strip(' \x00')
                      for ch in channels]
-            exclude = _find_exclude_idx(ch_names, exclude)
-            sel = list()
+            edf_info['units'] = list()
             for i, unit in enumerate(units):
+                if i in exclude:
+                    continue
                 if unit[:2] == 'uV':
-                    units[i] = 1e-6
+                    edf_info['units'].append(1e-6)
                 else:
-                    units[i] = 1
-                sel.append(i)
-            units = np.array(units, float)
+                    edf_info['units'].append(1)
+            edf_info['units'] = np.array(edf_info['units'], float)
 
             ch_names = [ch_names[idx] for idx in sel]
             physical_min = np.fromfile(fid, np.float64, len(channels))
@@ -762,8 +764,7 @@ def _read_gdf_header(fname, exclude):
                 meas_date=meas_date,
                 meas_id=meas_id, n_records=n_records, n_samps=n_samps,
                 nchan=nchan, subject_info=patient, physical_max=physical_max,
-                physical_min=physical_min, record_length=record_length,
-                units=units)
+                physical_min=physical_min, record_length=record_length)
 
             fid.seek(32 * edf_info['nchan'], 1)  # reserved
             assert fid.tell() == header_nbytes
@@ -902,6 +903,7 @@ def _read_gdf_header(fname, exclude):
             ch_names = [fid.read(16).decode().strip(' \x00')
                         for ch in channels]
             exclude = _find_exclude_idx(ch_names, exclude)
+            sel = np.setdiff1d(np.arange(len(ch_names)), exclude)
 
             fid.seek(80 * len(channels), 1)  # reserved space
             fid.seek(6 * len(channels), 1)  # phys_dim, obsolete
@@ -914,23 +916,24 @@ def _read_gdf_header(fname, exclude):
             """  # noqa
             units = np.fromfile(fid, np.uint16, len(channels)).tolist()
             unitcodes = np.array(units[:])
-            sel = list()
+            edf_info['units'] = list()
             for i, unit in enumerate(units):
+                if i in exclude:
+                    continue
                 if unit == 4275:  # microvolts
-                    units[i] = 1e-6
+                    edf_info['units'].append(1e-6)
                 elif unit == 4274:  # millivolts
-                    units[i] = 1e-3
+                    edf_info['units'].append(1e-3)
                 elif unit == 512:  # dimensionless
-                    units[i] = 1
+                    edf_info['units'].append(1)
                 elif unit == 0:
-                    units[i] = 1  # unrecognized
+                    edf_info['units'].append(1)  # unrecognized
                 else:
                     warn('Unsupported physical dimension for channel %d '
                          '(assuming dimensionless). Please contact the '
                          'MNE-Python developers for support.' % i)
-                    units[i] = 1
-                sel.append(i)
-            units = np.array(units, float)
+                    edf_info['units'].append(1)
+            edf_info['units'] = np.array(edf_info['units'], float)
 
             ch_names = [ch_names[idx] for idx in sel]
             physical_min = np.fromfile(fid, np.float64, len(channels))
@@ -989,7 +992,7 @@ def _read_gdf_header(fname, exclude):
                 meas_id=meas_id, n_records=n_records, n_samps=n_samps,
                 nchan=nchan, notch=notch, subject_info=patient,
                 physical_max=physical_max, physical_min=physical_min,
-                record_length=record_length, ref=ref, units=units)
+                record_length=record_length, ref=ref)
 
             # EVENT TABLE
             # -----------------------------------------------------------------

mne/io/edf/tests/test_gdf.py

@@ -116,4 +116,15 @@ def test_one_channel_gdf():
     assert 150.0 == ecg.info['sfreq']
 
 
+@testing.requires_testing_data
+def test_gdf_exclude_channels():
+    """Test reading GDF data with excluded channels."""
+    raw = read_raw_gdf(gdf1_path + '.gdf', exclude=('FP1', 'O1'))
+    assert 'FP1' not in raw.ch_names
+    assert 'O1' not in raw.ch_names
+    raw = read_raw_gdf(gdf2_path + '.gdf', exclude=('Fp1', 'O1'))
+    assert 'Fp1' not in raw.ch_names
+    assert 'O1' not in raw.ch_names
+
+
 run_tests_if_main()