Rainfall is poorly modeled by general circulation models (GCMs) and requires appropriate downscaling for local-scale hydrological impact studies. Such downscaling methods should be robust and accurate (to handle, e.g., extreme events and uncertainties), but the noncontinuous and highly nonlinear nature of rainfall makes this task particularly challenging. This paper brings together and extends state-of-the-art methods into an integrated and robust probabilistic methodology to downscale local daily rainfall series from an ensemble of climate simulations. The downscaling is based on generalized linear models (GLMs) that relate monthly GCM-scale atmospheric variables to local-scale daily rainfall series. A cross-validation step ensures that the fitted models are correctly conditioned by the climate variables, and a statistical procedure is proposed to test whether the statistical relationships identified for the reference period also hold in a future perturbed climate (i.e., to test the stationarity assumption). Additionally, we propose a strategy to downweigh poorly performing GCM-GLM couples. The methodology is assessed at 27 locations covering Switzerland and is shown to perform well in reproducing historical rainfall statistics including extremes and interannual variability. Furthermore, the projections are consistent with the simulations of physically based dynamical models. Using an original visualization method based on heat maps, we show that although the downscaling models were fitted at each of the 27 sites independently, their projections follow a spatially coherent pattern and that regions exhibiting different climate change impacts can be identified.