In this paper, we consider the near-critical dimer model in the setup of isoradial superpositions with Temperleyan boundary conditions. We show that the centered height function converges as the mesh size tends to zero to a limiting field which agrees with the (electromagnetically tilted) sine-Gordon model, whose derivative correlations are described by Grassmann variables (or equivalently determinants involving a massive Dirac operator). This answers a longstanding question in the field. A crucial part of the work is to develop a notion of discrete massive holomorphic functions and the tools to study such functions, in particular finding an exact discrete form of the massive Cauchy--Riemann equations, which is satisfied by the inverse Kasteleyn matrix. In comparison with previous studies, a key novelty of this part of our work is that the mass is not only allowed to be non-constant but can be complex-valued.