I have 2 kHz-sampled signal that I would like to downsample to 256 Hz
That's a resampling ratio of $\frac{256}{2000}=\frac{16}{125}$; not overly nice, but not too tragic to implement with a rational resampler, either.
If I understand well standard approaches to downsampling (and software implementing them) just plainly perform low pass and then selecting some bins from this new signal.
Well, this is not just downsampling, because 256 is not a factor of 2000; instead you'd theoretically need to upsample by 16, have an anti-image filter, an anti-alias filter, before downsampling by 125. Luckily, polyphase resampling has a method of doing that is mathematically identical to up- and then downconversion, but works at the lowest rate. Anyway, implementation details. At your rate, the data going through your signal processing chain is really negligible, so a naive implementation should do.
The problem is that when one does so, artifacts in the highres signal get propagated due to filtering and one sees them after downsampling in the time intervals that were "clean" before.
You will need to think about what "clean" really is, because you have a temporally well-limited problem, but that, conversely implies that any remedy that you apply to that region exclusively has wideband effects – effects on parts of the spectrum that were "clean" before.
is there an alternative way to perform "ecological" downsampling? E.g. if I just compute a moving average before selecting every n-th bin, would it make sense? (this way I'd get 2000 / 8 = 250 Hz signal, not 256 Hz, but it's ok for me)
That moving average is just a FIR filter, and not a very conveniently shaped one. Sure, that's easier to compute than a properly designed low-pass filter with arbitrary coefficients (instead of just "1" times the average length), but at a 2 kHz sampling rate, the mediocrest of 32bit microcontrollers could apply a pretty substantial FIR filter without breaking a sweat in real time. Don't think there's really a need to optimize this, if this is something running on something as powerful as a smartphone or laptop computer.
P.S. My data is MEG recordings on which I perform offline data analysis. Though I guess this should not matter here much.
If you're doing offline processing at so incredibly little data, then there's really no reason to avoid doing "good filters".
is there some open-source Python (or at least C/C++) package that performs "traditional" downsampling but taking care of artifact intervals?
You need to fix the intervals before you process the signal; that's probably not really hard to do. For example, you could replace these segments with zeros and let a low-pass filter run over the result to get rid of the windowing effects in spectrum. In fact, your resampling will necessitate filtering, anyway, so this might be "free".
How, however, you fix that artifact depends on the artifact and about the purpose of your signal processing. There can't be a "one size fits all" solution. So, no, there's no library that can do what you want - before you've identified what you need to do.
The resampling itself can be done by basically any DSP library I can think of: it's a rational resampler, with an interpolation of 16 and a decimation of 125.
(I think if you could ask another question that describes your artifacts, what "MEG" is, and what the purpose of your processing is, then we might be able to help you identify a method to solve the artifact problem. What you've asked here is interesting on its own, but its not really bringing you closer to a solution.)
