Seismic hazard assessment in active fault zones can benefit of strain rate measurements derived from geodetic data. Producing a continuous strain rate map from discrete data is an inverse problem traditionally tackled with standard interpolation schemes. Most algorithms require user-defined regression parameters that determine the smoothness of the recovered velocity field, and the amplitude of its spatial derivatives. This may lead to biases in the strain rates estimation which could eventually impact studies on earthquake hazard. Here we propose a transdimensional Bayesian method to estimate surface strain rates from GNSS velocities. We parameterize the velocity field with a variable number of Delaunay triangles, and use a reversible jump Monte-Carlo Markov Chain algorithm to sample the probability distribution of surface velocities and spatial derivatives. The solution is a complete probability distribution function for each component of the strain rate field. We conduct synthetic tests and compare our approach to a standard b-spline interpolation scheme. Our method is more resilient to data errors and uneven data distribution, while providing uncertainties associated with recovered velocities and strain rates. We apply our method to the Southwestern US, an extensively studied and monitored area and infer probabilistic strain rates along the main fault systems, including the San Andreas one, from the inversion of interseismic GNSS velocities. Our approach provide a full description of the strain rate tensor for zones where strain rates are highly contrasted, with no need to manually tune user-defined parameters. We recover sharp velocity gradients, without systematic artifacts.