Repository: mne-tools/mne-python

Issue:
Support STIM channels in EDF export
I found some issues with our recently added EDF export that we might want to fix before the next release:

1. The error if `EDFlib-Python` is not installed is misleading, because it says `RuntimeError: For exporting to EDF to work, the EDFlib module is needed, but it could not be imported`. However, although the import name is `EDFlib`, the PyPI name is actually `EDFlib-Python` and `EDFlib` is a different PyPI package (which requires Cython etc.).
2. `EDFlib-Python` is super slow compared to `pyEDFlib` (which uses the C implementation). I tested exporting with one of our BDF files from our lab (425MB), and it took `EDFlib-Python` 1 minute and 52 seconds to export. For comparison, `pyEDFlib` took only 6 seconds.
3. Analog stim channels are not exported. This is essential in BDF files, because Biosemi amplifiers do not have discrete triggers and instead always write an analog stim channel. There is no technical reason to not export stim channels in addition to EEG (other than carefully setting the scaling, unit, and physical min/max). I've fixed this in my EDF/BDF export implementation in MNELAB recently (today) in https://github.com/cbrnr/mnelab/pull/230.

The first issue is probably easy to fix. The second issue is pretty serious, we can't expect that users wait almost 2 minutes for something that can be done in 6 seconds. I'd consider switching to `pyEDFlib` instead (they have binary wheels and AFAIR the restriction to use pure Python applies only to our project but not external dependencies). The third issue requires some work, but this is important for BDF files.

Discussion:
welcome[bot]: Hello! 👋 Thanks for opening your first pull request here! ❤️ We will try to get back to you soon. 🚴

hofaflo: Conda+mamba checks fail since edfio is not on conda-forge yet, I'll work on that.
larsoner: > edfio is not on conda-forge yet, I'll work on that.

Let me know if you need help, I did this for `pysnirf2` just yesterday so it's fresh in my mind. But the short version is something like:

1. clone `staged-recipes`
2. `pip install grayskull` (or conda install it)
3. cd to `staged-recipes/recipes`
4. `grayskull pypi edfio`
5. Open PR to `staged-recipes`
hofaflo: Thank you for the simplified instructions!
Review is pending at conda-forge/staged-recipes#24511
Teuniz: 
>     * Two sections of `test_rawarray_edf` are removed, because
>       
>       * the patient birth date is formatted as "dd-MMM-yyyy" in EDF+, so a date like 20-JAN-1700 is actually valid
>       * the export always uses EDF+C, which allows to store dates outside the range valid in "original" Edf (1985-01-01 to 2084-12-31) formatted as "dd-MMM-yyyy" in the "local recording identification" field

Pay attention. Allowing a measurement date outside the range 1985-01-01 to 2084-12-31 breaks backwards compatibility with EDF!
In order to protect this backwards compatibility, EDFbrowser will refuse to open an EDF+ file that has a startdate like, for example, 1700-01-01 or 2090-01-01.
(Ofcourse, this limitation in EDFbrowser will be removed before the year 2085 arrives.)

hofaflo: Thank you for the heads up @Teuniz! edfio currently stores non-EDF-compatible startdates according to EDF+ specs in the "local recording identification" field and sets the EDF startdate field to "01.01.85". On read, it prefers the EDF+ startdate to the EDF startdate, but of course that only makes sense if no other software is used. Since EDFbrowser does not open files with such a mismatch, we will update edfio to refuse writing them :)
Teuniz: > Since EDFbrowser does not open files with such a mismatch, we will update edfio to refuse writing them :)

Thanks! In an "ideal" world where all EDF software has been updated to EDF+, this shouldn't be necessary,
but ofcourse this is not the case.
The goal is to avoid confusion when two different viewers (one EDF, the other EDF+) show different startdates.
Thanks again for your consideration. It's much appreciated.

hofaflo: [Conda installs](https://anaconda.org/conda-forge/edfio) should work now, but I'll mark this as draft until the-siesta-group/edfio#5 is merged and released.
sappelhoff: note: once this is merged, we'll have to make adjustments in mne-bids:

https://github.com/sappelhoff/mne-bids/blob/781fc577dc0da0419f7eb5c40d537b51834cb802/doc/install.rst?plain=1#L27
cbrnr: Just to be clear, even EDF+ mandates that the date range is between 1985–2084 (https://www.edfplus.info/specs/edfplus.html#additionalspecs). So regardless of EDFbrowser not opening it, this is non-standard and should therefore not be supported.

> (Ofcourse, this limitation in EDFbrowser will be removed before the year 2085 arrives.)

The spec also says: "After 2084, yy must be 'yy' and only item 4 of this paragraph defines the date."

And this actually made my day 😄!
cbrnr: Another comment: we don't actually require all dependencies to be on conda-forge, so installing via pip would have also worked. But I guess it doesn't hurt to have it available over there as well.
cbrnr: There's a problem: edfio requires Python >= 3.9, but MNE currently still supports Python 3.8.
hoechenberger: Mind you our installers are often lagging behind by 2 release or so, because of the complicated mix of dependencies. Otherwise I'm all for cutting off things that are holding us back. 
hofaflo: Is there already a timeline for dropping Python 3.8 support? If not – or if the plan is not to do that anytime soon – it should be relatively simple to support 3.8 in edfio. Let me know if that would be your preferred solution!


A few more things are fixed now:
- The dictionary mapping integers indicating sex in mne to letters used in EDF uses the mapping given [here](https://mne.tools/stable/generated/mne.Info.html#mne.Info)
- Empty subfields in patient/recording identification are substituted with `"X"` to enusre compatibility with EDF+
- Tolerances for comparing signal values and sample times to their values before exporting have been decreased to make passing the respective assertions non-trivial.
- Slicing signals to their original length has been removed where no padding occurs in the first place.

cbrnr: I'd prefer that we drop support for Python 3.8. This should be no problem, and we don't have to go full in and also drop 3.9. @drammock @larsoner OK?
larsoner: We try to follow the scientific Python community standards for this and I think that means it's okay to drop 3.8 according to a quick look at https://numpy.org/neps/nep-0029-deprecation_policy.html#support-table

Let's push 1.6 out the door Monday or Tuesday (it's overdue) then drop 3.8 it while bumping dependency versions to 2-years-old. In the meantime you can ignore the "old" 3.8 CI
cbrnr: > There's a problem: edfio requires Python >= 3.9, but MNE currently still supports Python 3.8.
> 
> This should be relatively easy to fix though by dropping support for Python 3.8. Actually, since we're already adopting SPECS, [SPEC 0](https://scientific-python.org/specs/spec-0000/) says we should drop 3.8 *and* 3.9 (and support only 3.10, 3.11, and 3.12).

SPEC0 is even more aggressive, but I'm fine with dropping only 3.8 after the 1.6 release.
cbrnr: One more thing before I'm happy: exporting stimulus channels are not supported, right? For example, [this test](https://github.com/mne-tools/mne-python/blob/main/mne/export/tests/test_export.py#L126) involves a `Raw` object with a stim channel, but it doesn't get exported. If this is intentional (I guess it is, but maybe we can take another look in another PR), I think we should at least emit a warning.
hofaflo: Good idea! Added a warning and removed the stim channel from all but one test to avoid checking for that warning in unrelated tests.
cbrnr: Great! Can you maybe remove [this comment](https://github.com/hofaflo/mne-python/blob/pr/hofaflo/12218/mne/export/tests/test_export.py#L181) (`# check channel types except for 'bio', which loses its type`), which is incorrect (`bio` stays `bio`)?

This test is a bit convoluted, feel free to clean up a bit if it makes sense. For example, there's actually a third export (the very last line), which I'm not sure what it actually tests (just a smoke test, but for what). And the test is called `test_double_export_edf()`, which is not true because of the third export...
Teuniz: > One more thing before I'm happy: exporting stimulus channels are not supported, right? For example, [this test](https://github.com/mne-tools/mne-python/blob/main/mne/export/tests/test_export.py#L126) involves a `Raw` object with a stim channel, but it doesn't get exported. If this is intentional (I guess it is, but maybe we can take another look in another PR), I think we should at least emit a warning.

Because I'm curious, what's different about a stim channel that it cannot be saved in EDF+?

cbrnr: > Because I'm curious, what's different about a stim channel that it cannot be saved in EDF+?

Good question, there's no reason why this cannot be done. In fact, this would be a perfectly valid and useful thing to do, but for some reason MNE-Python has excluded stim channels. We did not want to mess with that in this PR, but we'll tackle it in another PR. I guess the main difference between "normal" (physiological) and stim channels is that the latter contain exact (integer) values, but this can be dealt with by setting the physical and digital minima/maxima appropriately.
hofaflo: Deconvoluted the tests and rebased onto main, ready for re-review :)
hoechenberger: I think we'll need a changelog entry :)
hofaflo: Done!
cbrnr: Perfect, thanks!
cbrnr: I agree with your comments, so let's wait until they are addressed.

BTW, will this land in 1.6.1 (backported), or do we wait until 1.7? I'd prefer backporting.
hofaflo: Done, thanks for the suggestions @larsoner!

Backporting to 1.6.1 would require support for Python3.8 in edfio, right?
cbrnr: > Backporting to 1.6.1 would require support for Python3.8 in edfio, right?

True, I forgot. So I guess this will have to wait for v1.7 then (about 3 months).
hoechenberger: We should have na "Dependencies" section in the changelog to mention changes of dependencies. @larsoner WDYT?

Patch:
diff --git a/doc/changes/devel.rst b/doc/changes/devel.rst
index b30353407..e3738f86b 100644
--- a/doc/changes/devel.rst
+++ b/doc/changes/devel.rst
@@ -23,7 +23,8 @@ Version 1.7.dev0 (development)
 
 Enhancements
 ~~~~~~~~~~~~
-- None yet
+- Speed up export to .edf in :func:`mne.export.export_raw` by using ``edfio`` instead of ``EDFlib-Python`` (:gh:`12218` by :newcontrib:`Florian Hofer`)
+
 
 Bugs
 ~~~~
diff --git a/doc/changes/names.inc b/doc/changes/names.inc
index da884792c..2ec8f2268 100644
--- a/doc/changes/names.inc
+++ b/doc/changes/names.inc
@@ -172,6 +172,8 @@
 
 .. _Felix Raimundo: https://github.com/gamazeps
 
+.. _Florian Hofer: https://github.com/hofaflo
+
 .. _Florin Pop: https://github.com/florin-pop
 
 .. _Frederik Weber: https://github.com/Frederik-D-Weber
diff --git a/environment.yml b/environment.yml
index 3491974ff..75c57d693 100644
--- a/environment.yml
+++ b/environment.yml
@@ -56,7 +56,7 @@ dependencies:
   - mne-qt-browser
   - pymatreader
   - eeglabio
-  - edflib-python
+  - edfio>=0.2.1
   - pybv
   - mamba
   - lazy_loader
diff --git a/mne/export/_edf.py b/mne/export/_edf.py
index 7097f7bd8..04590f042 100644
--- a/mne/export/_edf.py
+++ b/mne/export/_edf.py
@@ -3,46 +3,15 @@
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
-from contextlib import contextmanager
+import datetime as dt
 
 import numpy as np
 
-from ..utils import _check_edflib_installed, warn
+from ..utils import _check_edfio_installed, warn
 
-_check_edflib_installed()
-from EDFlib.edfwriter import EDFwriter  # noqa: E402
-
-
-def _try_to_set_value(header, key, value, channel_index=None):
-    """Set key/value pairs in EDF header."""
-    # all EDFLib set functions are set<X>
-    # for example "setPatientName()"
-    func_name = f"set{key}"
-    func = getattr(header, func_name)
-
-    # some setter functions are indexed by channels
-    if channel_index is None:
-        return_val = func(value)
-    else:
-        return_val = func(channel_index, value)
-
-    # a nonzero return value indicates an error
-    if return_val != 0:
-        raise RuntimeError(
-            f"Setting {key} with {value} " f"returned an error value " f"{return_val}."
-        )
-
-
-@contextmanager
-def _auto_close(fid):
-    # try to close the handle no matter what
-    try:
-        yield fid
-    finally:
-        try:
-            fid.close()
-        except Exception:
-            pass  # we did our best
+_check_edfio_installed()
+from edfio import Edf, EdfAnnotation, EdfSignal, Patient, Recording  # noqa: E402
+from edfio._utils import round_float_to_8_characters  # noqa: E402
 
 
 def _export_raw(fname, raw, physical_range, add_ch_type):
@@ -51,9 +20,6 @@ def _export_raw(fname, raw, physical_range, add_ch_type):
     TODO: if in future the Info object supports transducer or
     technician information, allow writing those here.
     """
-    # scale to save data in EDF
-    phys_dims = "uV"
-
     # get EEG-related data in uV
     units = dict(
         eeg="uV", ecog="uV", seeg="uV", eog="uV", ecg="uV", emg="uV", bio="uV", dbs="uV"
@@ -61,7 +27,6 @@ def _export_raw(fname, raw, physical_range, add_ch_type):
 
     digital_min = -32767
     digital_max = 32767
-    file_type = EDFwriter.EDFLIB_FILETYPE_EDFPLUS
 
     # load data first
     raw.load_data()
@@ -73,6 +38,7 @@ def _export_raw(fname, raw, physical_range, add_ch_type):
         stim_index = np.argwhere(np.array(orig_ch_types) == "stim")
         stim_index = np.atleast_1d(stim_index.squeeze()).tolist()
         drop_chs.extend([raw.ch_names[idx] for idx in stim_index])
+        warn(f"Exporting STIM channels is not supported, dropping indices {stim_index}")
 
     # Add warning if any channel types are not voltage based.
     # Users are expected to only export data that is voltage based,
@@ -97,9 +63,11 @@ def _export_raw(fname, raw, physical_range, add_ch_type):
 
     ch_names = [ch for ch in raw.ch_names if ch not in drop_chs]
     ch_types = np.array(raw.get_channel_types(picks=ch_names))
-    n_channels = len(ch_names)
     n_times = raw.n_times
 
+    # get the entire dataset in uV
+    data = raw.get_data(units=units, picks=ch_names)
+
     # Sampling frequency in EDF only supports integers, so to allow for
     # float sampling rates from Raw, we adjust the output sampling rate
     # for all channels and the data record duration.
@@ -107,10 +75,20 @@ def _export_raw(fname, raw, physical_range, add_ch_type):
     if float(sfreq).is_integer():
         out_sfreq = int(sfreq)
         data_record_duration = None
+        # make non-integer second durations work
+        if (pad_width := int(np.ceil(n_times / sfreq) * sfreq - n_times)) > 0:
+            warn(
+                f"EDF format requires equal-length data blocks, "
+                f"so {pad_width / sfreq} seconds of "
+                "zeros were appended to all channels when writing the "
+                "final block."
+            )
+            data = np.pad(data, (0, int(pad_width)))
     else:
-        out_sfreq = np.floor(sfreq).astype(int)
-        data_record_duration = int(np.around(out_sfreq / sfreq, decimals=6) * 1e6)
-
+        data_record_duration = round_float_to_8_characters(
+            np.floor(sfreq) / sfreq, round
+        )
+        out_sfreq = np.floor(sfreq) / data_record_duration
         warn(
             f"Data has a non-integer sampling rate of {sfreq}; writing to "
             "EDF format may cause a small change to sample times."
@@ -122,16 +100,13 @@ def _export_raw(fname, raw, physical_range, add_ch_type):
     linefreq = raw.info["line_freq"]
     filter_str_info = f"HP:{highpass}Hz LP:{lowpass}Hz N:{linefreq}Hz"
 
-    # get the entire dataset in uV
-    data = raw.get_data(units=units, picks=ch_names)
-
     if physical_range == "auto":
         # get max and min for each channel type data
         ch_types_phys_max = dict()
         ch_types_phys_min = dict()
 
         for _type in np.unique(ch_types):
-            _picks = np.nonzero(ch_types == _type)[0]
+            _picks = [n for n, t in zip(ch_names, ch_types) if t == _type]
             _data = raw.get_data(units=units, picks=_picks)
             ch_types_phys_max[_type] = _data.max()
             ch_types_phys_min[_type] = _data.min()
@@ -157,178 +132,106 @@ def _export_raw(fname, raw, physical_range, add_ch_type):
                 f"The minimum μV of the data {data.min()} is "
                 f"less than the physical min passed in {pmin}.",
             )
+        data = np.clip(data, pmin, pmax)
+    signals = []
+    for idx, ch in enumerate(ch_names):
+        ch_type = ch_types[idx]
+        signal_label = f"{ch_type.upper()} {ch}" if add_ch_type else ch
+        if len(signal_label) > 16:
+            raise RuntimeError(
+                f"Signal label for {ch} ({ch_type}) is "
+                f"longer than 16 characters, which is not "
+                f"supported in EDF. Please shorten the "
+                f"channel name before exporting to EDF."
+            )
 
-    # create instance of EDF Writer
-    with _auto_close(EDFwriter(fname, file_type, n_channels)) as hdl:
-        # set channel data
-        for idx, ch in enumerate(ch_names):
-            ch_type = ch_types[idx]
-            signal_label = f"{ch_type.upper()} {ch}" if add_ch_type else ch
-            if len(signal_label) > 16:
-                raise RuntimeError(
-                    f"Signal label for {ch} ({ch_type}) is "
-                    f"longer than 16 characters, which is not "
-                    f"supported in EDF. Please shorten the "
-                    f"channel name before exporting to EDF."
-                )
-
-            if physical_range == "auto":
-                # take the channel type minimum and maximum
-                pmin = ch_types_phys_min[ch_type]
-                pmax = ch_types_phys_max[ch_type]
-            for key, val in [
-                ("PhysicalMaximum", pmax),
-                ("PhysicalMinimum", pmin),
-                ("DigitalMaximum", digital_max),
-                ("DigitalMinimum", digital_min),
-                ("PhysicalDimension", phys_dims),
-                ("SampleFrequency", out_sfreq),
-                ("SignalLabel", signal_label),
-                ("PreFilter", filter_str_info),
-            ]:
-                _try_to_set_value(hdl, key, val, channel_index=idx)
-
-        # set patient info
-        subj_info = raw.info.get("subject_info")
-        if subj_info is not None:
-            # get the full name of subject if available
-            first_name = subj_info.get("first_name", "")
-            middle_name = subj_info.get("middle_name", "")
-            last_name = subj_info.get("last_name", "")
-            name = " ".join(filter(None, [first_name, middle_name, last_name]))
-
-            birthday = subj_info.get("birthday")
-            hand = subj_info.get("hand")
-            weight = subj_info.get("weight")
-            height = subj_info.get("height")
-            sex = subj_info.get("sex")
-
-            additional_patient_info = []
-            for key, value in [("height", height), ("weight", weight), ("hand", hand)]:
-                if value:
-                    additional_patient_info.append(f"{key}={value}")
-            if len(additional_patient_info) == 0:
-                additional_patient_info = None
-            else:
-                additional_patient_info = " ".join(additional_patient_info)
-
-            if birthday is not None:
-                if hdl.setPatientBirthDate(birthday[0], birthday[1], birthday[2]) != 0:
-                    raise RuntimeError(
-                        f"Setting patient birth date to {birthday} "
-                        f"returned an error"
-                    )
-            for key, val in [
-                ("PatientCode", subj_info.get("his_id", "")),
-                ("PatientName", name),
-                ("PatientGender", sex),
-                ("AdditionalPatientInfo", additional_patient_info),
-            ]:
-                # EDFwriter compares integer encodings of sex and will
-                # raise a TypeError if value is None as returned by
-                # subj_info.get(key) if key is missing.
-                if val is not None:
-                    _try_to_set_value(hdl, key, val)
-
-        # set measurement date
-        meas_date = raw.info["meas_date"]
-        if meas_date:
-            subsecond = int(meas_date.microsecond / 100)
-            if (
-                hdl.setStartDateTime(
-                    year=meas_date.year,
-                    month=meas_date.month,
-                    day=meas_date.day,
-                    hour=meas_date.hour,
-                    minute=meas_date.minute,
-                    second=meas_date.second,
-                    subsecond=subsecond,
-                )
-                != 0
-            ):
-                raise RuntimeError(
-                    f"Setting start date time {meas_date} " f"returned an error"
-                )
-
-        device_info = raw.info.get("device_info")
-        if device_info is not None:
-            device_type = device_info.get("type")
-            _try_to_set_value(hdl, "Equipment", device_type)
-
-        # set data record duration
-        if data_record_duration is not None:
-            _try_to_set_value(hdl, "DataRecordDuration", data_record_duration)
-
-        # compute number of data records to loop over
-        n_blocks = np.ceil(n_times / out_sfreq).astype(int)
-
-        # increase the number of annotation signals if necessary
-        annots = raw.annotations
-        if annots is not None:
-            n_annotations = len(raw.annotations)
-            n_annot_chans = int(n_annotations / n_blocks) + 1
-            if n_annot_chans > 1:
-                hdl.setNumberOfAnnotationSignals(n_annot_chans)
-
-        # Write each data record sequentially
-        for idx in range(n_blocks):
-            end_samp = (idx + 1) * out_sfreq
-            if end_samp > n_times:
-                end_samp = n_times
-            start_samp = idx * out_sfreq
-
-            # then for each datarecord write each channel
-            for jdx in range(n_channels):
-                # create a buffer with sampling rate
-                buf = np.zeros(out_sfreq, np.float64, "C")
+        if physical_range == "auto":
+            # take the channel type minimum and maximum
+            pmin = ch_types_phys_min[ch_type]
+            pmax = ch_types_phys_max[ch_type]
+
+        signals.append(
+            EdfSignal(
+                data[idx],
+                out_sfreq,
+                label=signal_label,
+                transducer_type="",
+                physical_dimension="uV",
+                physical_range=(pmin, pmax),
+                digital_range=(digital_min, digital_max),
+                prefiltering=filter_str_info,
+            )
+        )
 
-                # get channel data for this block
-                ch_data = data[jdx, start_samp:end_samp]
+    # set patient info
+    subj_info = raw.info.get("subject_info")
+    if subj_info is not None:
+        # get the full name of subject if available
+        first_name = subj_info.get("first_name", "")
+        middle_name = subj_info.get("middle_name", "")
+        last_name = subj_info.get("last_name", "")
+        name = "_".join(filter(None, [first_name, middle_name, last_name]))
+
+        birthday = subj_info.get("birthday")
+        if birthday is not None:
+            birthday = dt.date(*birthday)
+        hand = subj_info.get("hand")
+        weight = subj_info.get("weight")
+        height = subj_info.get("height")
+        sex = subj_info.get("sex")
+
+        additional_patient_info = []
+        for key, value in [("height", height), ("weight", weight), ("hand", hand)]:
+            if value:
+                additional_patient_info.append(f"{key}={value}")
+
+        patient = Patient(
+            code=subj_info.get("his_id") or "X",
+            sex={0: "X", 1: "M", 2: "F", None: "X"}[sex],
+            birthdate=birthday,
+            name=name or "X",
+            additional=additional_patient_info,
+        )
+    else:
+        patient = None
 
-                # assign channel data to the buffer and write to EDF
-                buf[: len(ch_data)] = ch_data
-                err = hdl.writeSamples(buf)
-                if err != 0:
-                    raise RuntimeError(
-                        f"writeSamples() for channel{ch_names[jdx]} "
-                        f"returned error: {err}"
-                    )
+    # set measurement date
+    if (meas_date := raw.info["meas_date"]) is not None:
+        startdate = dt.date(meas_date.year, meas_date.month, meas_date.day)
+        starttime = dt.time(
+            meas_date.hour, meas_date.minute, meas_date.second, meas_date.microsecond
+        )
+    else:
+        startdate = None
+        starttime = None
 
-            # there was an incomplete datarecord
-            if len(ch_data) != len(buf):
-                warn(
-                    f"EDF format requires equal-length data blocks, "
-                    f"so {(len(buf) - len(ch_data)) / sfreq} seconds of "
-                    "zeros were appended to all channels when writing the "
-                    "final block."
+    device_info = raw.info.get("device_info")
+    if device_info is not None:
+        device_type = device_info.get("type") or "X"
+        recording = Recording(startdate=startdate, equipment_code=device_type)
+    else:
+        recording = Recording(startdate=startdate)
+
+    annotations = []
+    for desc, onset, duration, ch_names in zip(
+        raw.annotations.description,
+        raw.annotations.onset,
+        raw.annotations.duration,
+        raw.annotations.ch_names,
+    ):
+        if ch_names:
+            for ch_name in ch_names:
+                annotations.append(
+                    EdfAnnotation(onset, duration, desc + f"@@{ch_name}")
                 )
-
-        # write annotations
-        if annots is not None:
-            for desc, onset, duration, ch_names in zip(
-                raw.annotations.description,
-                raw.annotations.onset,
-                raw.annotations.duration,
-                raw.annotations.ch_names,
-            ):
-                # annotations are written in terms of 100 microseconds
-                onset = onset * 10000
-                duration = duration * 10000
-                if ch_names:
-                    for ch_name in ch_names:
-                        if (
-                            hdl.writeAnnotation(onset, duration, desc + f"@@{ch_name}")
-                            != 0
-                        ):
-                            raise RuntimeError(
-                                f"writeAnnotation() returned an error "
-                                f"trying to write {desc}@@{ch_name} at {onset} "
-                                f"for {duration} seconds."
-                            )
-                else:
-                    if hdl.writeAnnotation(onset, duration, desc) != 0:
-                        raise RuntimeError(
-                            f"writeAnnotation() returned an error "
-                            f"trying to write {desc} at {onset} "
-                            f"for {duration} seconds."
-                        )
+        else:
+            annotations.append(EdfAnnotation(onset, duration, desc))
+
+    Edf(
+        signals=signals,
+        patient=patient,
+        recording=recording,
+        starttime=starttime,
+        data_record_duration=data_record_duration,
+        annotations=annotations,
+    ).write(fname)
diff --git a/mne/utils/__init__.pyi b/mne/utils/__init__.pyi
index 42694921f..3e4d1292e 100644
--- a/mne/utils/__init__.pyi
+++ b/mne/utils/__init__.pyi
@@ -32,7 +32,7 @@ __all__ = [
     "_check_depth",
     "_check_dict_keys",
     "_check_dt",
-    "_check_edflib_installed",
+    "_check_edfio_installed",
     "_check_eeglabio_installed",
     "_check_event_id",
     "_check_fname",
@@ -230,7 +230,7 @@ from .check import (
     _check_compensation_grade,
     _check_depth,
     _check_dict_keys,
-    _check_edflib_installed,
+    _check_edfio_installed,
     _check_eeglabio_installed,
     _check_event_id,
     _check_fname,
diff --git a/mne/utils/check.py b/mne/utils/check.py
index eb8e14de2..8c2bc5f91 100644
--- a/mne/utils/check.py
+++ b/mne/utils/check.py
@@ -11,11 +11,9 @@ import re
 from builtins import input  # no-op here but facilitates testing
 from difflib import get_close_matches
 from importlib import import_module
-from importlib.metadata import version
 from pathlib import Path
 
 import numpy as np
-from packaging.version import parse
 
 from ..defaults import HEAD_SIZE_DEFAULT, _handle_default
 from ..fixes import _compare_version, _median_complex
@@ -368,7 +366,6 @@ def _soft_import(name, purpose, strict=True):
     # Mapping import namespaces to their pypi package name
     pip_name = dict(
         sklearn="scikit-learn",
-        EDFlib="EDFlib-Python",
         mne_bids="mne-bids",
         mne_nirs="mne-nirs",
         mne_features="mne-features",
@@ -411,21 +408,9 @@ def _check_eeglabio_installed(strict=True):
     return _soft_import("eeglabio", "exporting to EEGLab", strict=strict)
 
 
-def _check_edflib_installed(strict=True):
+def _check_edfio_installed(strict=True):
     """Aux function."""
-    out = _soft_import("EDFlib", "exporting to EDF", strict=strict)
-    if out:
-        # EDFlib-Python 1.0.7 is not compatible with NumPy 2.0
-        # https://gitlab.com/Teuniz/EDFlib-Python/-/issues/10
-        ver = version("EDFlib-Python")
-        if parse(ver) <= parse("1.0.7") and parse(np.__version__).major >= 2:
-            if strict:  # pragma: no cover
-                raise RuntimeError(
-                    f"EDFlib version={ver} is not compatible with NumPy 2.0, consider "
-                    "upgrading EDFlib-Python"
-                )
-            out = False
-    return out
+    return _soft_import("edfio", "exporting to EDF", strict=strict)
 
 
 def _check_pybv_installed(strict=True):
diff --git a/mne/utils/docs.py b/mne/utils/docs.py
index d8eb668ae..015a5ff7d 100644
--- a/mne/utils/docs.py
+++ b/mne/utils/docs.py
@@ -1463,7 +1463,7 @@ docdict[
 Supported formats:
     - BrainVision (``.vhdr``, ``.vmrk``, ``.eeg``, uses `pybv <https://github.com/bids-standard/pybv>`_)
     - EEGLAB (``.set``, uses :mod:`eeglabio`)
-    - EDF (``.edf``, uses `EDFlib-Python <https://gitlab.com/Teuniz/EDFlib-Python>`_)
+    - EDF (``.edf``, uses `edfio <https://github.com/the-siesta-group/edfio>`_)
 """  # noqa: E501
 
 docdict[
diff --git a/pyproject.toml b/pyproject.toml
index 1ff5d9313..9c21ef711 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -122,7 +122,7 @@ test_extra = [
     "nbclient",
     "sphinx-gallery",
     "eeglabio",
-    "EDFlib-Python",
+    "edfio>=0.2.1",
     "pybv",
     "imageio>=2.6.1",
     "imageio-ffmpeg>=0.4.1",


Test patch:
diff --git a/mne/export/tests/test_export.py b/mne/export/tests/test_export.py
index 67bd417bb..4e86c3bb6 100644
--- a/mne/export/tests/test_export.py
+++ b/mne/export/tests/test_export.py
@@ -6,7 +6,6 @@
 
 from contextlib import nullcontext
 from datetime import datetime, timezone
-from os import remove
 from pathlib import Path
 
 import numpy as np
@@ -33,7 +32,7 @@ from mne.io import (
 )
 from mne.tests.test_epochs import _get_data
 from mne.utils import (
-    _check_edflib_installed,
+    _check_edfio_installed,
     _record_warnings,
     _resource_path,
     object_diff,
@@ -120,17 +119,11 @@ def test_export_raw_eeglab(tmp_path):
         raw.export(temp_fname, overwrite=True)
 
 
-@pytest.mark.skipif(
-    not _check_edflib_installed(strict=False), reason="edflib-python not installed"
-)
-def test_double_export_edf(tmp_path):
-    """Test exporting an EDF file multiple times."""
-    rng = np.random.RandomState(123456)
-    format = "edf"
+def _create_raw_for_edf_tests(stim_channel_index=None):
+    rng = np.random.RandomState(12345)
     ch_types = [
         "eeg",
         "eeg",
-        "stim",
         "ecog",
         "ecog",
         "seeg",
@@ -140,12 +133,27 @@ def test_double_export_edf(tmp_path):
         "dbs",
         "bio",
     ]
-    info = create_info(len(ch_types), sfreq=1000, ch_types=ch_types)
-    info = info.set_meas_date("2023-09-04 14:53:09.000")
-    data = rng.random(size=(len(ch_types), 1000)) * 1e-5
+    if stim_channel_index is not None:
+        ch_types.insert(stim_channel_index, "stim")
+    ch_names = np.arange(len(ch_types)).astype(str).tolist()
+    info = create_info(ch_names, sfreq=1000, ch_types=ch_types)
+    data = rng.random(size=(len(ch_names), 2000)) * 1e-5
+    return RawArray(data, info)
+
+
+edfio_mark = pytest.mark.skipif(
+    not _check_edfio_installed(strict=False), reason="edfio not installed"
+)
+
+
+@edfio_mark()
+def test_double_export_edf(tmp_path):
+    """Test exporting an EDF file multiple times."""
+    raw = _create_raw_for_edf_tests(stim_channel_index=2)
+    raw.info.set_meas_date("2023-09-04 14:53:09.000")
 
     # include subject info and measurement date
-    info["subject_info"] = dict(
+    raw.info["subject_info"] = dict(
         his_id="12345",
         first_name="mne",
         last_name="python",
@@ -155,15 +163,14 @@ def test_double_export_edf(tmp_path):
         height=1.75,
         hand=3,
     )
-    raw = RawArray(data, info)
 
     # export once
-    temp_fname = tmp_path / f"test.{format}"
-    raw.export(temp_fname, add_ch_type=True)
+    temp_fname = tmp_path / "test.edf"
+    with pytest.warns(RuntimeWarning, match="Exporting STIM channels"):
+        raw.export(temp_fname, add_ch_type=True)
     raw_read = read_raw_edf(temp_fname, infer_types=True, preload=True)
 
     # export again
-    raw_read.load_data()
     raw_read.export(temp_fname, add_ch_type=True, overwrite=True)
     raw_read = read_raw_edf(temp_fname, infer_types=True, preload=True)
 
@@ -171,53 +178,22 @@ def test_double_export_edf(tmp_path):
     raw.drop_channels("2")
 
     assert raw.ch_names == raw_read.ch_names
-    # only compare the original length, since extra zeros are appended
-    orig_raw_len = len(raw)
-    assert_array_almost_equal(
-        raw.get_data(), raw_read.get_data()[:, :orig_raw_len], decimal=4
-    )
-    assert_allclose(raw.times, raw_read.times[:orig_raw_len], rtol=0, atol=1e-5)
+    assert_array_almost_equal(raw.get_data(), raw_read.get_data(), decimal=10)
+    assert_array_equal(raw.times, raw_read.times)
 
     # check info
     for key in set(raw.info) - {"chs"}:
         assert raw.info[key] == raw_read.info[key]
 
-    # check channel types except for 'bio', which loses its type
     orig_ch_types = raw.get_channel_types()
     read_ch_types = raw_read.get_channel_types()
     assert_array_equal(orig_ch_types, read_ch_types)
 
-    # check handling of missing subject metadata
-    del info["subject_info"]["sex"]
-    raw_2 = RawArray(data, info)
-    raw_2.export(temp_fname, add_ch_type=True, overwrite=True)
 
-
-@pytest.mark.skipif(
-    not _check_edflib_installed(strict=False), reason="edflib-python not installed"
-)
+@edfio_mark()
 def test_export_edf_annotations(tmp_path):
     """Test that exporting EDF preserves annotations."""
-    rng = np.random.RandomState(123456)
-    format = "edf"
-    ch_types = [
-        "eeg",
-        "eeg",
-        "stim",
-        "ecog",
-        "ecog",
-        "seeg",
-        "eog",
-        "ecg",
-        "emg",
-        "dbs",
-        "bio",
-    ]
-    ch_names = np.arange(len(ch_types)).astype(str).tolist()
-    info = create_info(ch_names, sfreq=1000, ch_types=ch_types)
-    data = rng.random(size=(len(ch_names), 2000)) * 1.0e-5
-    raw = RawArray(data, info)
-
+    raw = _create_raw_for_edf_tests()
     annotations = Annotations(
         onset=[0.01, 0.05, 0.90, 1.05],
         duration=[0, 1, 0, 0],
@@ -227,7 +203,7 @@ def test_export_edf_annotations(tmp_path):
     raw.set_annotations(annotations)
 
     # export
-    temp_fname = tmp_path / f"test.{format}"
+    temp_fname = tmp_path / "test.edf"
     raw.export(temp_fname)
 
     # read in the file
@@ -238,24 +214,19 @@ def test_export_edf_annotations(tmp_path):
     assert_array_equal(raw.annotations.ch_names, raw_read.annotations.ch_names)
 
 
-@pytest.mark.skipif(
-    not _check_edflib_installed(strict=False), reason="edflib-python not installed"
-)
+@edfio_mark()
 def test_rawarray_edf(tmp_path):
     """Test saving a Raw array with integer sfreq to EDF."""
-    rng = np.random.RandomState(12345)
-    format = "edf"
-    ch_types = ["eeg", "eeg", "stim", "ecog", "seeg", "eog", "ecg", "emg", "dbs", "bio"]
-    ch_names = np.arange(len(ch_types)).astype(str).tolist()
-    info = create_info(ch_names, sfreq=1000, ch_types=ch_types)
-    data = rng.random(size=(len(ch_names), 1000)) * 1e-5
+    raw = _create_raw_for_edf_tests()
 
     # include subject info and measurement date
-    subject_info = dict(
-        first_name="mne", last_name="python", birthday=(1992, 1, 20), sex=1, hand=3
+    raw.info["subject_info"] = dict(
+        first_name="mne",
+        last_name="python",
+        birthday=(1992, 1, 20),
+        sex=1,
+        hand=3,
     )
-    info["subject_info"] = subject_info
-    raw = RawArray(data, info)
     time_now = datetime.now()
     meas_date = datetime(
         year=time_now.year,
@@ -267,125 +238,104 @@ def test_rawarray_edf(tmp_path):
         tzinfo=timezone.utc,
     )
     raw.set_meas_date(meas_date)
-    temp_fname = tmp_path / f"test.{format}"
+    temp_fname = tmp_path / "test.edf"
 
     raw.export(temp_fname, add_ch_type=True)
     raw_read = read_raw_edf(temp_fname, infer_types=True, preload=True)
 
-    # stim channel should be dropped
-    raw.drop_channels("2")
-
     assert raw.ch_names == raw_read.ch_names
-    # only compare the original length, since extra zeros are appended
-    orig_raw_len = len(raw)
-    assert_array_almost_equal(
-        raw.get_data(), raw_read.get_data()[:, :orig_raw_len], decimal=4
-    )
-    assert_allclose(raw.times, raw_read.times[:orig_raw_len], rtol=0, atol=1e-5)
+    assert_array_almost_equal(raw.get_data(), raw_read.get_data(), decimal=10)
+    assert_array_equal(raw.times, raw_read.times)
 
-    # check channel types except for 'bio', which loses its type
     orig_ch_types = raw.get_channel_types()
     read_ch_types = raw_read.get_channel_types()
     assert_array_equal(orig_ch_types, read_ch_types)
     assert raw.info["meas_date"] == raw_read.info["meas_date"]
 
-    # channel name can't be longer than 16 characters with the type added
-    raw_bad = raw.copy()
-    raw_bad.rename_channels({"1": "abcdefghijklmnopqrstuvwxyz"})
-    with pytest.raises(RuntimeError, match="Signal label"), pytest.warns(
-        RuntimeWarning, match="Data has a non-integer"
-    ):
-        raw_bad.export(temp_fname, overwrite=True)
 
-    # include bad birthday that is non-EDF compliant
-    bad_info = info.copy()
-    bad_info["subject_info"]["birthday"] = (1700, 1, 20)
-    raw = RawArray(data, bad_info)
-    with pytest.raises(RuntimeError, match="Setting patient birth date"):
-        raw.export(temp_fname, overwrite=True)
+@edfio_mark()
+def test_edf_export_warns_on_non_voltage_channels(tmp_path):
+    """Test saving a Raw array containing a non-voltage channel."""
+    temp_fname = tmp_path / "test.edf"
 
-    # include bad measurement date that is non-EDF compliant
-    raw = RawArray(data, info)
-    meas_date = datetime(year=1984, month=1, day=1, tzinfo=timezone.utc)
-    raw.set_meas_date(meas_date)
-    with pytest.raises(RuntimeError, match="Setting start date time"):
-        raw.export(temp_fname, overwrite=True)
-
-    # test that warning is raised if there are non-voltage based channels
-    raw = RawArray(data, info)
+    raw = _create_raw_for_edf_tests()
     raw.set_channel_types({"9": "hbr"}, on_unit_change="ignore")
     with pytest.warns(RuntimeWarning, match="Non-voltage channels"):
         raw.export(temp_fname, overwrite=True)
 
     # data should match up to the non-accepted channel
     raw_read = read_raw_edf(temp_fname, preload=True)
-    orig_raw_len = len(raw)
-    assert_array_almost_equal(
-        raw.get_data()[:-1, :], raw_read.get_data()[:, :orig_raw_len], decimal=4
-    )
-    assert_allclose(raw.times, raw_read.times[:orig_raw_len], rtol=0, atol=1e-5)
-
-    # the data should still match though
-    raw_read = read_raw_edf(temp_fname, preload=True)
-    raw.drop_channels("2")
     assert raw.ch_names == raw_read.ch_names
-    orig_raw_len = len(raw)
-    assert_array_almost_equal(
-        raw.get_data(), raw_read.get_data()[:, :orig_raw_len], decimal=4
-    )
-    assert_allclose(raw.times, raw_read.times[:orig_raw_len], rtol=0, atol=1e-5)
+    assert_array_almost_equal(raw.get_data()[:-1], raw_read.get_data()[:-1], decimal=10)
+    assert_array_equal(raw.times, raw_read.times)
+
+
+@edfio_mark()
+def test_channel_label_too_long_for_edf_raises_error(tmp_path):
+    """Test trying to save an EDF where a channel label is longer than 16 characters."""
+    raw = _create_raw_for_edf_tests()
+    raw.rename_channels({"1": "abcdefghijklmnopqrstuvwxyz"})
+    with pytest.raises(RuntimeError, match="Signal label"):
+        raw.export(tmp_path / "test.edf")
 
 
-@pytest.mark.skipif(
-    not _check_edflib_installed(strict=False), reason="edflib-python not installed"
+@edfio_mark()
+def test_measurement_date_outside_range_valid_for_edf(tmp_path):
+    """Test trying to save an EDF with a measurement date before 1985-01-01."""
+    raw = _create_raw_for_edf_tests()
+    raw.set_meas_date(datetime(year=1984, month=1, day=1, tzinfo=timezone.utc))
+    with pytest.raises(ValueError, match="EDF only allows dates from 1985 to 2084"):
+        raw.export(tmp_path / "test.edf", overwrite=True)
+
+
+@pytest.mark.parametrize(
+    ("physical_range", "exceeded_bound"),
+    [
+        ((-1e6, 0), "maximum"),
+        ((0, 1e6), "minimum"),
+    ],
 )
+@edfio_mark()
+def test_export_edf_signal_clipping(tmp_path, physical_range, exceeded_bound):
+    """Test if exporting data exceeding physical min/max clips and emits a warning."""
+    raw = read_raw_fif(fname_raw)
+    raw.pick(picks=["eeg", "ecog", "seeg"]).load_data()
+    temp_fname = tmp_path / "test.edf"
+    with pytest.warns(RuntimeWarning, match=f"The {exceeded_bound}"):
+        raw.export(temp_fname, physical_range=physical_range)
+    raw_read = read_raw_edf(temp_fname, preload=True)
+    assert raw_read.get_data().min() >= physical_range[0]
+    assert raw_read.get_data().max() <= physical_range[1]
+
+
+@edfio_mark()
 @pytest.mark.parametrize(
-    ["dataset", "format"],
+    ("input_path", "warning_msg"),
     [
-        ["test", "edf"],
-        pytest.param("misc", "edf", marks=[pytest.mark.slowtest, misc._pytest_mark()]),
+        (fname_raw, "Data has a non-integer"),
+        pytest.param(
+            misc_path / "ecog" / "sample_ecog_ieeg.fif",
+            "EDF format requires",
+            marks=[pytest.mark.slowtest, misc._pytest_mark()],
+        ),
     ],
 )
-def test_export_raw_edf(tmp_path, dataset, format):
+def test_export_raw_edf(tmp_path, input_path, warning_msg):
     """Test saving a Raw instance to EDF format."""
-    if dataset == "test":
-        raw = read_raw_fif(fname_raw)
-    elif dataset == "misc":
-        fname = misc_path / "ecog" / "sample_ecog_ieeg.fif"
-        raw = read_raw_fif(fname)
+    raw = read_raw_fif(input_path)
 
     # only test with EEG channels
     raw.pick(picks=["eeg", "ecog", "seeg"]).load_data()
-    orig_ch_names = raw.ch_names
-    temp_fname = tmp_path / f"test.{format}"
-
-    # test runtime errors
-    with pytest.warns() as record:
-        raw.export(temp_fname, physical_range=(-1e6, 0))
-    if dataset == "test":
-        assert any("Data has a non-integer" in str(rec.message) for rec in record)
-    assert any("The maximum" in str(rec.message) for rec in record)
-    remove(temp_fname)
-
-    with pytest.warns() as record:
-        raw.export(temp_fname, physical_range=(0, 1e6))
-    if dataset == "test":
-        assert any("Data has a non-integer" in str(rec.message) for rec in record)
-    assert any("The minimum" in str(rec.message) for rec in record)
-    remove(temp_fname)
-
-    if dataset == "test":
-        with pytest.warns(RuntimeWarning, match="Data has a non-integer"):
-            raw.export(temp_fname)
-    elif dataset == "misc":
-        with pytest.warns(RuntimeWarning, match="EDF format requires"):
-            raw.export(temp_fname)
+    temp_fname = tmp_path / "test.edf"
+
+    with pytest.warns(RuntimeWarning, match=warning_msg):
+        raw.export(temp_fname)
 
     if "epoc" in raw.ch_names:
         raw.drop_channels(["epoc"])
 
     raw_read = read_raw_edf(temp_fname, preload=True)
-    assert orig_ch_names == raw_read.ch_names
+    assert raw.ch_names == raw_read.ch_names
     # only compare the original length, since extra zeros are appended
     orig_raw_len = len(raw)
 
@@ -395,7 +345,7 @@ def test_export_raw_edf(tmp_path, dataset, format):
     # will result in a resolution of 0.09 uV. This resolution
     # though is acceptable for most EEG manufacturers.
     assert_array_almost_equal(
-        raw.get_data(), raw_read.get_data()[:, :orig_raw_len], decimal=4
+        raw.get_data(), raw_read.get_data()[:, :orig_raw_len], decimal=8
     )
 
     # Due to the data record duration limitations of EDF files, one
@@ -407,6 +357,27 @@ def test_export_raw_edf(tmp_path, dataset, format):
     assert_allclose(raw.times, raw_read.times[:orig_raw_len], rtol=0, atol=1e-5)
 
 
+@edfio_mark()
+def test_export_raw_edf_does_not_fail_on_empty_header_fields(tmp_path):
+    """Test writing a Raw instance with empty header fields to EDF."""
+    rng = np.random.RandomState(123456)
+
+    ch_types = ["eeg"]
+    info = create_info(len(ch_types), sfreq=1000, ch_types=ch_types)
+    info["subject_info"] = {
+        "his_id": "",
+        "first_name": "",
+        "middle_name": "",
+        "last_name": "",
+    }
+    info["device_info"] = {"type": "123"}
+
+    data = rng.random(size=(len(ch_types), 1000)) * 1e-5
+    raw = RawArray(data, info)
+
+    raw.export(tmp_path / "test.edf", add_ch_type=True)
+
+
 @pytest.mark.xfail(reason="eeglabio (usage?) bugs that should be fixed")
 @pytest.mark.parametrize("preload", (True, False))
 def test_export_epochs_eeglab(tmp_path, preload):