  We apply for the first time Gaussian Process Regression (GPR) as a foreground
removal technique in the context of single-dish, low redshift HI intensity
mapping, and present an open-source Python toolkit for doing so. We use MeerKAT
and SKA1-MID-like simulations of 21cm foregrounds (including polarisation
leakage), HI cosmological signal and instrumental noise. We find that it is
possible to use GPR as a foreground removal technique in this context, and that
it is better suited in some cases to recover the HI power spectrum than
Principal Component Analysis (PCA), especially on small scales. GPR is
especially good at recovering the radial power spectrum, outperforming PCA when
considering the full bandwidth of our data. Both methods are worse at
recovering the transverse power spectrum, since they rely on frequency-only
covariance information. When halving our data along frequency, we find that GPR
performs better in the low frequency range, where foregrounds are brighter. It
performs worse than PCA when frequency channels are missing, to emulate RFI
flagging. We conclude that GPR is an excellent foreground removal option for
the case of single-dish, low redshift HI intensity mapping in the absence of
missing frequency channels. Our Python toolkit gpr4im and the data used in this
analysis are publicly available on GitHub.
