The Monte Carlo method is used to simulate diffusion of interacting lattice gases on heterogeneous surfaces. The fluctuation and Kubo-Green methods are utilized for determining the tracer and chemical surface diffusion coefficients. Simulations are carried out on a square array of 64 x 64 lattice sites applying periodic boundary conditions. Surface heterogeneity is introduced in terms of the bivariate trap model with random topography. Simulations are performed for a trap concentration theta(trap) = 0.2 and various values of the trap binding energies epsilon(2), epsilon(1) and the nearest neighbor interaction energy phi(NN). For phi(NN) < 0 (repulsion) the system exhibits c(2 x 2) ordering at half coverage and low temperatures. Due to the low compressibility of the well-ordered c(2 x 2) lattice gas phase, the chemical diffusion coefficient is large under these circumstances. The effects of surface heterogeneity, epsilon(1) not equal epsilon(2), are largely pronounced at low total coverages theta and low temperatures T where most of the adatoms are trapped by the deep traps. At higher coverages the effects are relatively small albeit significant. The existence of deep traps disturbs c(2 x 2) ordering at theta = 0.5 and low temperatures. It is shown that the breakdown of the order substantially affects surface diffusion.