Surface water waves are considered propagating over highly variable non-smooth topographies. For this three-dimensional problem a Dirichlet-to-Neumann (DtN) operator is constructed reducing the numerical modelling and evolution to the two-dimensional free surface. The corresponding discrete Fourier integral operator is defined through a matrix decomposition. The topographic component of the decomposition requires special care, and a Galerkin method is provided accordingly. One-dimensional numerical simulations, along the free surface, validate the DtN formulation in the presence of a large-amplitude rapidly varying topography. An alternative conformal-mapping-based method is used for benchmarking. A two-dimensional simulation in the presence of a Luneburg lens (a particular submerged mound) illustrates the accurate performance of the three-dimensional DtN operator.