Radio weak lensing, while a highly promising complementary probe to optical weak lensing, will require incredible precision in the measurement of galaxy shape parameters. In this paper, we extend the Bayesian Inference for Radio Observations model fitting approach to measure galaxy shapes directly from visibility data of radio continuum surveys, instead of from image data. We apply a Hamiltonian Monte Carlo (HMC) technique for sampling the posterior, which is more efficient than the standard Monte Carlo Markov Chain method when dealing with a large dimensional parameter space. Adopting the exponential profile for galaxy model fitting allows us to analytically calculate the likelihood gradient required by HMC, allowing a faster and more accurate sampling. The method is tested on SKA1-MID simulated observations at 1.4 GHz of a field containing up to 1000 star-forming galaxies. It is also applied to a simulated observation of the weak lensing precursor survey SuperCLASS. In both cases we obtain reliable measurements of the galaxies' ellipticity and size for all sources with SNR $\ge 10$, and we also find relationships between the convergence properties of the HMC technique and some source parameters. Direct shape measurement in the visibility domain achieves high accuracy at the expected source number densities of the current and next SKA precursor continuum surveys. The proposed method can be easily extended for the fitting of other galaxy and scientific parameters, as well as simultaneously marginalising over systematic and instrumental effects.